src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/PairCollision.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2011 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>
  33 Foam::scalar Foam::PairCollision<CloudType>::cosPhiMinFlatWall = 1 - SMALL;
  34
  35 template<class CloudType>
  36 Foam::scalar Foam::PairCollision<CloudType>::flatWallDuplicateExclusion =
  37     sqrt(3*SMALL);
  38
  39
  40 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  41
  42 template<class CloudType>
  43 void Foam::PairCollision<CloudType>::preInteraction()
  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>
  58 void Foam::PairCollision<CloudType>::parcelInteraction()
  59 {
  60     PstreamBuffers pBufs(Pstream::nonBlocking);
  61
  62     il_.sendReferredData(this->owner().cellOccupancy(), pBufs);
  63
  64     realRealInteraction();
  65
  66     il_.receiveReferredData(pBufs);
  67
  68     realReferredInteraction();
  69 }
  70
  71
  72 template<class CloudType>
  73 void Foam::PairCollision<CloudType>::realRealInteraction()
  74 {
  75     // Direct interaction list (dil)
  76     const labelListList& dil = il_.dil();
  77
  78     typename CloudType::parcelType* pA_ptr = NULL;
  79     typename CloudType::parcelType* pB_ptr = NULL;
  80
  81     List<DynamicList<typename CloudType::parcelType*> >& cellOccupancy =
  82         this->owner().cellOccupancy();
  83
  84     forAll(dil, realCellI)
  85     {
  86         // Loop over all Parcels in cell A (a)
  87         forAll(cellOccupancy[realCellI], a)
  88         {
  89             pA_ptr = cellOccupancy[realCellI][a];
  90
  91             forAll(dil[realCellI], interactingCells)
  92             {
  93                 List<typename CloudType::parcelType*> cellBParcels =
  94                     cellOccupancy[dil[realCellI][interactingCells]];
  95
  96                 // Loop over all Parcels in cell B (b)
  97                 forAll(cellBParcels, b)
  98                 {
  99                     pB_ptr = cellBParcels[b];
 100
 101                     evaluatePair(*pA_ptr, *pB_ptr);
 102                 }
 103             }
 104
 105             // Loop over the other Parcels in cell A (aO)
 106             forAll(cellOccupancy[realCellI], aO)
 107             {
 108                 pB_ptr = cellOccupancy[realCellI][aO];
 109
 110                 // Do not double-evaluate, compare pointers, arbitrary
 111                 // order
 112                 if (pB_ptr > pA_ptr)
 113                 {
 114                     evaluatePair(*pA_ptr, *pB_ptr);
 115                 }
 116             }
 117         }
 118     }
 119 }
 120
 121
 122 template<class CloudType>
 123 void Foam::PairCollision<CloudType>::realReferredInteraction()
 124 {
 125     // Referred interaction list (ril)
 126     const labelListList& ril = il_.ril();
 127
 128     List<IDLList<typename CloudType::parcelType> >& referredParticles =
 129     il_.referredParticles();
 130
 131     List<DynamicList<typename CloudType::parcelType*> >& cellOccupancy =
 132         this->owner().cellOccupancy();
 133
 134     // Loop over all referred cells
 135     forAll(ril, refCellI)
 136     {
 137         IDLList<typename CloudType::parcelType>& refCellRefParticles =
 138             referredParticles[refCellI];
 139
 140         const labelList& realCells = ril[refCellI];
 141
 142         // Loop over all referred parcels in the referred cell
 143
 144         forAllIter
 145         (
 146             typename IDLList<typename CloudType::parcelType>,
 147             refCellRefParticles,
 148             referredParcel
 149         )
 150         {
 151             // Loop over all real cells in that the referred cell is
 152             // to supply interactions to
 153
 154             forAll(realCells, realCellI)
 155             {
 156                 List<typename CloudType::parcelType*> realCellParcels =
 157                     cellOccupancy[realCells[realCellI]];
 158
 159                 forAll(realCellParcels, realParcelI)
 160                 {
 161                     evaluatePair
 162                     (
 163                         *realCellParcels[realParcelI],
 164                         referredParcel()
 165                     );
 166                 }
 167             }
 168         }
 169     }
 170 }
 171
 172
 173 template<class CloudType>
 174 void Foam::PairCollision<CloudType>::wallInteraction()
 175 {
 176     const polyMesh& mesh = this->owner().mesh();
 177
 178     const labelListList& dil = il_.dil();
 179
 180     const labelListList directWallFaces = il_.dwfil();
 181
 182     const labelList& patchID = mesh.boundaryMesh().patchID();
 183
 184     const volVectorField& U = mesh.lookupObject<volVectorField>(il_.UName());
 185
 186     List<DynamicList<typename CloudType::parcelType*> >& cellOccupancy =
 187         this->owner().cellOccupancy();
 188
 189     // Storage for the wall interaction sites
 190     DynamicList<point> flatSitePoints;
 191     DynamicList<scalar> flatSiteExclusionDistancesSqr;
 192     DynamicList<WallSiteData<vector> > flatSiteData;
 193     DynamicList<point> otherSitePoints;
 194     DynamicList<scalar> otherSiteDistances;
 195     DynamicList<WallSiteData<vector> > otherSiteData;
 196     DynamicList<point> sharpSitePoints;
 197     DynamicList<scalar> sharpSiteExclusionDistancesSqr;
 198     DynamicList<WallSiteData<vector> > sharpSiteData;
 199
 200     forAll(dil, realCellI)
 201     {
 202         // The real wall faces in range of this real cell
 203         const labelList& realWallFaces = directWallFaces[realCellI];
 204
 205         // Loop over all Parcels in cell
 206         forAll(cellOccupancy[realCellI], cellParticleI)
 207         {
 208             flatSitePoints.clear();
 209             flatSiteExclusionDistancesSqr.clear();
 210             flatSiteData.clear();
 211             otherSitePoints.clear();
 212             otherSiteDistances.clear();
 213             otherSiteData.clear();
 214             sharpSitePoints.clear();
 215             sharpSiteExclusionDistancesSqr.clear();
 216             sharpSiteData.clear();
 217
 218             typename CloudType::parcelType& p =
 219                 *cellOccupancy[realCellI][cellParticleI];
 220
 221             const point& pos = p.position();
 222
 223             scalar r = wallModel_->pREff(p);
 224
 225             // real wallFace interactions
 226
 227             forAll(realWallFaces, realWallFaceI)
 228             {
 229                 label realFaceI = realWallFaces[realWallFaceI];
 230
 231                 pointHit nearest = mesh.faces()[realFaceI].nearestPoint
 232                 (
 233                     pos,
 234                     mesh.points()
 235                 );
 236
 237                 if (nearest.distance() < r)
 238                 {
 239                     vector normal = mesh.faceAreas()[realFaceI];
 240
 241                     normal /= mag(normal);
 242
 243                     const vector& nearPt = nearest.rawPoint();
 244
 245                     vector pW = nearPt - pos;
 246
 247                     scalar normalAlignment = normal & pW/mag(pW);
 248
 249                     // Find the patchIndex and wallData for WallSiteData object
 250                     label patchI = patchID[realFaceI - mesh.nInternalFaces()];
 251
 252                     label patchFaceI =
 253                         realFaceI - mesh.boundaryMesh()[patchI].start();
 254
 255                     WallSiteData<vector> wSD
 256                     (
 257                         patchI,
 258                         U.boundaryField()[patchI][patchFaceI]
 259                     );
 260
 261                     bool particleHit = false;
 262                     if (normalAlignment > cosPhiMinFlatWall)
 263                     {
 264                         // Guard against a flat interaction being
 265                         // present on the boundary of two or more
 266                         // faces, which would create duplicate contact
 267                         // points. Duplicates are discarded.
 268                         if
 269                         (
 270                             !duplicatePointInList
 271                             (
 272                                 flatSitePoints,
 273                                 nearPt,
 274                                 sqr(r*flatWallDuplicateExclusion)
 275                             )
 276                         )
 277                         {
 278                             flatSitePoints.append(nearPt);
 279
 280                             flatSiteExclusionDistancesSqr.append
 281                             (
 282                                 sqr(r) - sqr(nearest.distance())
 283                             );
 284
 285                             flatSiteData.append(wSD);
 286
 287                             particleHit = true;
 288                         }
 289                     }
 290                     else
 291                     {
 292                         otherSitePoints.append(nearPt);
 293
 294                         otherSiteDistances.append(nearest.distance());
 295
 296                         otherSiteData.append(wSD);
 297
 298                         particleHit = true;
 299                     }
 300
 301                     if (particleHit)
 302                     {
 303                         this->owner().functions().postFace(p, realFaceI);
 304                         this->owner().functions().postPatch
 305                         (
 306                             p,
 307                             patchI,
 308                             patchFaceI
 309                         );
 310                     }
 311                 }
 312             }
 313
 314             // referred wallFace interactions
 315
 316             // The labels of referred wall faces in range of this real cell
 317             const labelList& cellRefWallFaces = il_.rwfilInverse()[realCellI];
 318
 319
 320
 321             forAll(cellRefWallFaces, rWFI)
 322             {
 323                 label refWallFaceI = cellRefWallFaces[rWFI];
 324
 325                 const referredWallFace& rwf =
 326                     il_.referredWallFaces()[refWallFaceI];
 327
 328                 const pointField& pts = rwf.points();
 329
 330                 pointHit nearest = rwf.nearestPoint(pos, pts);
 331
 332                 if (nearest.distance() < r)
 333                 {
 334                     vector normal = rwf.normal(pts);
 335
 336                     normal /= mag(normal);
 337
 338                     const vector& nearPt = nearest.rawPoint();
 339
 340                     vector pW = nearPt - pos;
 341
 342                     scalar normalAlignment = normal & pW/mag(pW);
 343
 344                     // Find the patchIndex and wallData for WallSiteData object
 345
 346                     WallSiteData<vector> wSD
 347                     (
 348                         rwf.patchIndex(),
 349                         il_.referredWallData()[refWallFaceI]
 350                     );
 351
 352                     bool particleHit = false;
 353                     if (normalAlignment > cosPhiMinFlatWall)
 354                     {
 355                         // Guard against a flat interaction being
 356                         // present on the boundary of two or more
 357                         // faces, which would create duplicate contact
 358                         // points. Duplicates are discarded.
 359                         if
 360                         (
 361                             !duplicatePointInList
 362                             (
 363                                 flatSitePoints,
 364                                 nearPt,
 365                                 sqr(r*flatWallDuplicateExclusion)
 366                             )
 367                         )
 368                         {
 369                             flatSitePoints.append(nearPt);
 370
 371                             flatSiteExclusionDistancesSqr.append
 372                             (
 373                                 sqr(r) - sqr(nearest.distance())
 374                             );
 375
 376                             flatSiteData.append(wSD);
 377
 378                             particleHit = true;
 379                         }
 380                     }
 381                     else
 382                     {
 383                         otherSitePoints.append(nearPt);
 384
 385                         otherSiteDistances.append(nearest.distance());
 386
 387                         otherSiteData.append(wSD);
 388
 389                         particleHit = true;
 390                     }
 391
 392                     if (particleHit)
 393                     {
 394                         // TODO: call cloud function objects for referred
 395                         //       wall particle interactions
 396                     }
 397                 }
 398             }
 399
 400             // All flat interaction sites found, now classify the
 401             // other sites as being in range of a flat interaction, or
 402             // a sharp interaction, being aware of not duplicating the
 403             // sharp interaction sites.
 404
 405             // The "other" sites need to evaluated in order of
 406             // ascending distance to their nearest point so that
 407             // grouping occurs around the closest in any group
 408
 409             labelList sortedOtherSiteIndices;
 410
 411             sortedOrder(otherSiteDistances, sortedOtherSiteIndices);
 412
 413             forAll(sortedOtherSiteIndices, siteI)
 414             {
 415                 label orderedIndex = sortedOtherSiteIndices[siteI];
 416
 417                 const point& otherPt = otherSitePoints[orderedIndex];
 418
 419                 if
 420                 (
 421                     !duplicatePointInList
 422                     (
 423                         flatSitePoints,
 424                         otherPt,
 425                         flatSiteExclusionDistancesSqr
 426                     )
 427                 )
 428                 {
 429                     // Not in range of a flat interaction, must be a
 430                     // sharp interaction.
 431
 432                     if
 433                     (
 434                         !duplicatePointInList
 435                         (
 436                             sharpSitePoints,
 437                             otherPt,
 438                             sharpSiteExclusionDistancesSqr
 439                         )
 440                     )
 441                     {
 442                         sharpSitePoints.append(otherPt);
 443
 444                         sharpSiteExclusionDistancesSqr.append
 445                         (
 446                             sqr(r) - sqr(otherSiteDistances[orderedIndex])
 447                         );
 448
 449                         sharpSiteData.append(otherSiteData[orderedIndex]);
 450                     }
 451                 }
 452             }
 453
 454             evaluateWall
 455             (
 456                 p,
 457                 flatSitePoints,
 458                 flatSiteData,
 459                 sharpSitePoints,
 460                 sharpSiteData
 461             );
 462         }
 463     }
 464 }
 465
 466
 467 template<class CloudType>
 468 bool Foam::PairCollision<CloudType>::duplicatePointInList
 469 (
 470     const DynamicList<point>& existingPoints,
 471     const point& pointToTest,
 472     scalar duplicateRangeSqr
 473 ) const
 474 {
 475     forAll(existingPoints, i)
 476     {
 477         if (magSqr(existingPoints[i] - pointToTest) < duplicateRangeSqr)
 478         {
 479             return true;
 480         }
 481     }
 482
 483     return false;
 484 }
 485
 486
 487 template<class CloudType>
 488 bool Foam::PairCollision<CloudType>::duplicatePointInList
 489 (
 490     const DynamicList<point>& existingPoints,
 491     const point& pointToTest,
 492     const scalarList& duplicateRangeSqr
 493 ) const
 494 {
 495     forAll(existingPoints, i)
 496     {
 497         if (magSqr(existingPoints[i] - pointToTest) < duplicateRangeSqr[i])
 498         {
 499             return true;
 500         }
 501     }
 502
 503     return false;
 504 }
 505
 506
 507 template<class CloudType>
 508 void Foam::PairCollision<CloudType>::postInteraction()
 509 {
 510     // Delete any collision records where no collision occurred this step
 511
 512     forAllIter(typename CloudType, this->owner(), iter)
 513     {
 514         typename CloudType::parcelType& p = iter();
 515
 516         p.collisionRecords().update();
 517     }
 518 }
 519
 520
 521 template<class CloudType>
 522 void Foam::PairCollision<CloudType>::evaluatePair
 523 (
 524     typename CloudType::parcelType& pA,
 525     typename CloudType::parcelType& pB
 526 ) const
 527 {
 528     pairModel_->evaluatePair(pA, pB);
 529 }
 530
 531
 532 template<class CloudType>
 533 void Foam::PairCollision<CloudType>::evaluateWall
 534 (
 535     typename CloudType::parcelType& p,
 536     const List<point>& flatSitePoints,
 537     const List<WallSiteData<vector> >& flatSiteData,
 538     const List<point>& sharpSitePoints,
 539     const List<WallSiteData<vector> >& sharpSiteData
 540 ) const
 541 {
 542     wallModel_->evaluateWall
 543     (
 544         p,
 545         flatSitePoints,
 546         flatSiteData,
 547         sharpSitePoints,
 548         sharpSiteData
 549     );
 550 }
 551
 552
 553 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 554
 555 template<class CloudType>
 556 Foam::PairCollision<CloudType>::PairCollision
 557 (
 558     const dictionary& dict,
 559     CloudType& owner
 560 )
 561 :
 562     CollisionModel<CloudType>(dict, owner, typeName),
 563     pairModel_
 564     (
 565         PairModel<CloudType>::New
 566         (
 567             this->coeffDict(),
 568             this->owner()
 569         )
 570     ),
 571     wallModel_
 572     (
 573         WallModel<CloudType>::New
 574         (
 575             this->coeffDict(),
 576             this->owner()
 577         )
 578     ),
 579     il_
 580     (
 581         owner.mesh(),
 582         readScalar(this->coeffDict().lookup("maxInteractionDistance")),
 583         Switch
 584         (
 585             this->coeffDict().lookupOrDefault
 586             (
 587                 "writeReferredParticleCloud",
 588                 false
 589             )
 590         ),
 591         this->coeffDict().lookupOrDefault("UName", word("U"))
 592     )
 593 {}
 594
 595
 596 template<class CloudType>
 597 Foam::PairCollision<CloudType>::PairCollision(PairCollision<CloudType>& cm)
 598 :
 599     CollisionModel<CloudType>(cm),
 600     pairModel_(NULL),
 601     wallModel_(NULL),
 602     il_(cm.owner().mesh())
 603 {
 604     notImplemented
 605     (
 606         "Foam::PairCollision<CloudType>::PairCollision"
 607         "("
 608             "PairCollision<CloudType>& cm"
 609         ")"
 610     );
 611 }
 612
 613
 614 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 615
 616 template<class CloudType>
 617 Foam::PairCollision<CloudType>::~PairCollision()
 618 {}
 619
 620
 621 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 622
 623 template<class CloudType>
 624 Foam::label Foam::PairCollision<CloudType>::nSubCycles() const
 625 {
 626     label nSubCycles = 1;
 627
 628     if (pairModel_->controlsTimestep())
 629     {
 630         label nPairSubCycles = returnReduce
 631         (
 632             pairModel_->nSubCycles(), maxOp<label>()
 633         );
 634
 635         nSubCycles = max(nSubCycles, nPairSubCycles);
 636     }
 637
 638     if (wallModel_->controlsTimestep())
 639     {
 640         label nWallSubCycles = returnReduce
 641         (
 642             wallModel_->nSubCycles(), maxOp<label>()
 643         );
 644
 645         nSubCycles = max(nSubCycles, nWallSubCycles);
 646     }
 647
 648     return nSubCycles;
 649 }
 650
 651
 652 template<class CloudType>
 653 bool Foam::PairCollision<CloudType>::controlsWallInteraction() const
 654 {
 655     return true;
 656 }
 657
 658
 659 template<class CloudType>
 660 void Foam::PairCollision<CloudType>::collide()
 661 {
 662     preInteraction();
 663
 664     parcelInteraction();
 665
 666     wallInteraction();
 667
 668     postInteraction();
 669 }
 670
 671
 672 // ************************************************************************* //