src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2009-2009 OpenCFD Ltd.
   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 the
  13     Free Software Foundation; either version 2 of the License, or (at your
  14     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, write to the Free Software Foundation,
  23     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
  24
  25 \*---------------------------------------------------------------------------*/
  26
  27 #include "DsmcCloud.H"
  28 #include "BinaryCollisionModel.H"
  29 #include "WallInteractionModel.H"
  30 #include "InflowBoundaryModel.H"
  31
  32 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  33
  34 template<class ParcelType>
  35 Foam::scalar Foam::DsmcCloud<ParcelType>::kb = 1.380650277e-23;
  36
  37
  38 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  39
  40 template<class ParcelType>
  41 void Foam::DsmcCloud<ParcelType>::buildConstProps()
  42 {
  43     Info<< nl << "Constructing constant properties for" << endl;
  44     constProps_.setSize(typeIdList_.size());
  45
  46     dictionary moleculeProperties
  47     (
  48         particleProperties_.subDict("moleculeProperties")
  49     );
  50
  51     forAll(typeIdList_, i)
  52     {
  53         const word& id(typeIdList_[i]);
  54
  55         Info<< "    " << id << endl;
  56
  57         const dictionary& molDict(moleculeProperties.subDict(id));
  58
  59         constProps_[i] =
  60         typename ParcelType::constantProperties::constantProperties(molDict);
  61     }
  62 }
  63
  64
  65 template<class ParcelType>
  66 void Foam::DsmcCloud<ParcelType>::buildCellOccupancy()
  67 {
  68     forAll(cellOccupancy_, cO)
  69     {
  70         cellOccupancy_[cO].clear();
  71     }
  72
  73     forAllIter(typename DsmcCloud<ParcelType>, *this, iter)
  74     {
  75         cellOccupancy_[iter().cell()].append(&iter());
  76     }
  77 }
  78
  79
  80 template<class ParcelType>
  81 void Foam::DsmcCloud<ParcelType>::initialise
  82 (
  83     const IOdictionary& dsmcInitialiseDict
  84 )
  85 {
  86     Info<< nl << "Initialising particles" << endl;
  87
  88     const scalar temperature
  89     (
  90         readScalar(dsmcInitialiseDict.lookup("temperature"))
  91     );
  92
  93     const vector velocity(dsmcInitialiseDict.lookup("velocity"));
  94
  95     const dictionary& numberDensitiesDict
  96     (
  97         dsmcInitialiseDict.subDict("numberDensities")
  98     );
  99
 100     List<word> molecules(numberDensitiesDict.toc());
 101
 102     Field<scalar> numberDensities(molecules.size());
 103
 104     forAll(molecules, i)
 105     {
 106         numberDensities[i] = readScalar
 107         (
 108             numberDensitiesDict.lookup(molecules[i])
 109         );
 110     }
 111
 112     numberDensities /= nParticle_;
 113
 114     forAll(mesh_.cells(), cell)
 115     {
 116         const vector& cC = mesh_.cellCentres()[cell];
 117         const labelList& cellFaces = mesh_.cells()[cell];
 118         const scalar cV = mesh_.cellVolumes()[cell];
 119
 120         label nTets = 0;
 121
 122         // Each face is split into nEdges (or nVertices) - 2 tets.
 123         forAll(cellFaces, face)
 124         {
 125             nTets += mesh_.faces()[cellFaces[face]].size() - 2;
 126         }
 127
 128         // Calculate the cumulative tet volumes circulating around the cell and
 129         // record the vertex labels of each.
 130         scalarList cTetVFracs(nTets, 0.0);
 131
 132         List<labelList> tetPtIs(nTets, labelList(3,-1));
 133
 134         // Keep track of which tet this is.
 135         label tet = 0;
 136
 137         forAll(cellFaces, face)
 138         {
 139             const labelList& facePoints = mesh_.faces()[cellFaces[face]];
 140
 141             label pointI = 1;
 142             while ((pointI + 1) < facePoints.size())
 143             {
 144
 145                 const vector& pA = mesh_.points()[facePoints[0]];
 146                 const vector& pB = mesh_.points()[facePoints[pointI]];
 147                 const vector& pC = mesh_.points()[facePoints[pointI + 1]];
 148
 149                 cTetVFracs[tet] =
 150                     mag(((pA - cC) ^ (pB - cC)) & (pC - cC))/(cV*6.0)
 151                   + cTetVFracs[max((tet - 1),0)];
 152
 153                 tetPtIs[tet][0] = facePoints[0];
 154                 tetPtIs[tet][1] = facePoints[pointI];
 155                 tetPtIs[tet][2] = facePoints[pointI + 1];
 156
 157                 pointI++;
 158                 tet++;
 159             }
 160         }
 161
 162         // Force the last volume fraction value to 1.0 to avoid any
 163         // rounding/non-flat face errors giving a value < 1.0
 164         cTetVFracs[nTets - 1] = 1.0;
 165
 166         forAll(molecules, i)
 167         {
 168             const word& moleculeName(molecules[i]);
 169
 170             label typeId(findIndex(typeIdList_, moleculeName));
 171
 172             if (typeId == -1)
 173             {
 174                 FatalErrorIn("Foam::DsmcCloud<ParcelType>::initialise")
 175                 << "typeId " << moleculeName << "not defined." << nl
 176                     << abort(FatalError);
 177             }
 178
 179             const typename ParcelType::constantProperties& cP =
 180                 constProps(typeId);
 181
 182             scalar numberDensity = numberDensities[i];
 183
 184             // Calculate the number of particles required
 185             scalar particlesRequired = numberDensity*mesh_.cellVolumes()[cell];
 186
 187             // Only integer numbers of particles can be inserted
 188             label nParticlesToInsert = label(particlesRequired);
 189
 190             // Add another particle with a probability proportional to the
 191             // remainder of taking the integer part of particlesRequired
 192             if ((particlesRequired - nParticlesToInsert) > rndGen_.scalar01())
 193             {
 194                 nParticlesToInsert++;
 195             }
 196
 197             for (label pI = 0; pI < nParticlesToInsert; pI++)
 198             {
 199                 // Choose a random point in a generic tetrahedron
 200
 201                 scalar s = rndGen_.scalar01();
 202                 scalar t = rndGen_.scalar01();
 203                 scalar u = rndGen_.scalar01();
 204
 205                 if (s + t > 1.0)
 206                 {
 207                     s = 1.0 - s;
 208                     t = 1.0 - t;
 209                 }
 210
 211                 if (t + u > 1.0)
 212                 {
 213                     scalar tmp = u;
 214                     u = 1.0 - s - t;
 215                     t = 1.0 - tmp;
 216                 }
 217                 else if (s + t + u > 1.0)
 218                 {
 219                     scalar tmp = u;
 220                     u = s + t + u - 1.0;
 221                     s = 1.0 - t - tmp;
 222                 }
 223
 224                 // Choose a tetrahedron to insert in, based on their relative
 225                 // volumes
 226                 scalar tetSelection = rndGen_.scalar01();
 227
 228                 // Selected tetrahedron
 229                 label sTet = -1;
 230
 231                 forAll(cTetVFracs, tet)
 232                 {
 233                     sTet = tet;
 234
 235                     if (cTetVFracs[tet] >= tetSelection)
 236                     {
 237                         break;
 238                     }
 239                 }
 240
 241                 vector p =
 242                     (1 - s - t - u)*cC
 243                   + s*mesh_.points()[tetPtIs[sTet][0]]
 244                   + t*mesh_.points()[tetPtIs[sTet][1]]
 245                   + u*mesh_.points()[tetPtIs[sTet][2]];
 246
 247                 vector U = equipartitionLinearVelocity
 248                 (
 249                     temperature,
 250                     cP.mass()
 251                 );
 252
 253                 scalar Ei = equipartitionInternalEnergy
 254                 (
 255                     temperature,
 256                     cP.internalDegreesOfFreedom()
 257                 );
 258
 259                 U += velocity;
 260
 261                 addNewParcel
 262                 (
 263                     p,
 264                     U,
 265                     Ei,
 266                     cell,
 267                     typeId
 268                 );
 269             }
 270         }
 271     }
 272
 273     // Initialise the sigmaTcRMax_ field to the product of the cross section of
 274     // the most abundant species and the most probable thermal speed (Bird,
 275     // p222-223)
 276
 277     label mostAbundantType(findMax(numberDensities));
 278
 279     const typename ParcelType::constantProperties& cP = constProps
 280     (
 281         mostAbundantType
 282     );
 283
 284     sigmaTcRMax_.internalField() = cP.sigmaT()*maxwellianMostProbableSpeed
 285     (
 286         temperature,
 287         cP.mass()
 288     );
 289
 290     sigmaTcRMax_.correctBoundaryConditions();
 291 }
 292
 293
 294 template<class ParcelType>
 295 void Foam::DsmcCloud<ParcelType>::collisions()
 296 {
 297     buildCellOccupancy();
 298
 299     // Temporary storage for subCells
 300     List<DynamicList<label> > subCells(8);
 301
 302     scalar deltaT = mesh().time().deltaTValue();
 303
 304     label collisionCandidates = 0;
 305
 306     label collisions = 0;
 307
 308     forAll(cellOccupancy_, celli)
 309     {
 310         const DynamicList<ParcelType*>& cellParcels(cellOccupancy_[celli]);
 311
 312         label nC(cellParcels.size());
 313
 314         if (nC > 1)
 315         {
 316
 317             // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 318             // Assign particles to one of 8 Cartesian subCells
 319
 320             // Clear temporary lists
 321             forAll(subCells, i)
 322             {
 323                 subCells[i].clear();
 324             }
 325
 326             // Inverse addressing specifying which subCell a parcel is in
 327             List<label> whichSubCell(cellParcels.size());
 328
 329             const point& cC = mesh_.cellCentres()[celli];
 330
 331             forAll(cellParcels, i)
 332             {
 333                 ParcelType* p = cellParcels[i];
 334
 335                 vector relPos = p->position() - cC;
 336
 337                 label subCell =
 338                     pos(relPos.x()) + 2*pos(relPos.y()) + 4*pos(relPos.z());
 339
 340                 subCells[subCell].append(i);
 341
 342                 whichSubCell[i] = subCell;
 343             }
 344
 345             // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 346
 347             scalar sigmaTcRMax = sigmaTcRMax_[celli];
 348
 349             scalar selectedPairs =
 350                 collisionSelectionRemainder_[celli]
 351               + 0.5*nC*(nC - 1)*nParticle_*sigmaTcRMax*deltaT
 352                /mesh_.cellVolumes()[celli];
 353
 354             label nCandidates(selectedPairs);
 355
 356             collisionSelectionRemainder_[celli] = selectedPairs - nCandidates;
 357
 358             collisionCandidates += nCandidates;
 359
 360             for (label c = 0; c < nCandidates; c++)
 361             {
 362                 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 363                 // subCell candidate selection procedure
 364
 365                 // Select the first collision candidate
 366                 label candidateP = rndGen_.integer(0, nC - 1);
 367
 368                 // Declare the second collision candidate
 369                 label candidateQ = -1;
 370
 371                 List<label> subCellPs = subCells[whichSubCell[candidateP]];
 372
 373                 label nSC = subCellPs.size();
 374
 375                 if (nSC > 1)
 376                 {
 377                     // If there are two or more particle in a subCell, choose
 378                     // another from the same cell.  If the same candidate is
 379                     // chosen, choose again.
 380
 381                     do
 382                     {
 383                         candidateQ = subCellPs[rndGen_.integer(0, nSC - 1)];
 384
 385                     } while(candidateP == candidateQ);
 386                 }
 387                 else
 388                 {
 389                     // Select a possible second collision candidate from the
 390                     // whole cell.  If the same candidate is chosen, choose
 391                     // again.
 392
 393                     do
 394                     {
 395                         candidateQ = rndGen_.integer(0, nC - 1);
 396
 397                     } while(candidateP == candidateQ);
 398                 }
 399
 400                 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 401                 // uniform candidate selection procedure
 402
 403                 // // Select the first collision candidate
 404                 // label candidateP = rndGen_.integer(0, nC-1);
 405
 406                 // // Select a possible second collision candidate
 407                 // label candidateQ = rndGen_.integer(0, nC-1);
 408
 409                 // // If the same candidate is chosen, choose again
 410                 // while(candidateP == candidateQ)
 411                 // {
 412                 //     candidateQ = rndGen_.integer(0, nC-1);
 413                 // }
 414
 415                 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 416
 417                 ParcelType* parcelP = cellParcels[candidateP];
 418                 ParcelType* parcelQ = cellParcels[candidateQ];
 419
 420                 label typeIdP = parcelP->typeId();
 421                 label typeIdQ = parcelQ->typeId();
 422
 423                 scalar sigmaTcR = binaryCollision().sigmaTcR
 424                 (
 425                     typeIdP,
 426                     typeIdQ,
 427                     parcelP->U(),
 428                     parcelQ->U()
 429                 );
 430
 431                 // Update the maximum value of sigmaTcR stored, but use the
 432                 // initial value in the acceptance-rejection criteria because
 433                 // the number of collision candidates selected was based on this
 434
 435                 if (sigmaTcR > sigmaTcRMax_[celli])
 436                 {
 437                     sigmaTcRMax_[celli] = sigmaTcR;
 438                 }
 439
 440                 if ((sigmaTcR/sigmaTcRMax) > rndGen_.scalar01())
 441                 {
 442                     binaryCollision().collide
 443                     (
 444                         typeIdP,
 445                         typeIdQ,
 446                         parcelP->U(),
 447                         parcelQ->U(),
 448                         parcelP->Ei(),
 449                         parcelQ->Ei()
 450                     );
 451
 452                     collisions++;
 453                 }
 454             }
 455         }
 456     }
 457
 458     reduce(collisions, sumOp<label>());
 459
 460     reduce(collisionCandidates, sumOp<label>());
 461
 462     if (collisionCandidates)
 463     {
 464         Info<< "    Collisions                      = "
 465             << collisions << nl
 466             << "    Acceptance rate                 = "
 467             << scalar(collisions)/scalar(collisionCandidates) << nl
 468             << endl;
 469     }
 470     else
 471     {
 472         Info<< "    No collisions" << endl;
 473     }
 474 }
 475
 476
 477 template<class ParcelType>
 478 void Foam::DsmcCloud<ParcelType>::resetFields()
 479 {
 480     q_ = dimensionedScalar("zero",  dimensionSet(1, 0, -3, 0, 0), 0.0);
 481
 482     fD_ = dimensionedVector
 483     (
 484         "zero",
 485         dimensionSet(1, -1, -2, 0, 0),
 486         vector::zero
 487     );
 488
 489     rhoN_ = dimensionedScalar("zero",  dimensionSet(0, -3, 0, 0, 0), VSMALL);
 490
 491     rhoM_ =  dimensionedScalar("zero",  dimensionSet(1, -3, 0, 0, 0), VSMALL);
 492
 493     dsmcRhoN_ = dimensionedScalar("zero",  dimensionSet(0, -3, 0, 0, 0), 0.0);
 494
 495     linearKE_ = dimensionedScalar("zero",  dimensionSet(1, -1, -2, 0, 0), 0.0);
 496
 497     internalE_ = dimensionedScalar("zero",  dimensionSet(1, -1, -2, 0, 0), 0.0);
 498
 499     iDof_ = dimensionedScalar("zero",  dimensionSet(0, -3, 0, 0, 0), VSMALL);
 500
 501     momentum_ = dimensionedVector
 502     (
 503         "zero",
 504         dimensionSet(1, -2, -1, 0, 0),
 505         vector::zero
 506     );
 507 }
 508
 509
 510 template<class ParcelType>
 511 void Foam::DsmcCloud<ParcelType>::calculateFields()
 512 {
 513     scalarField& rhoN = rhoN_.internalField();
 514
 515     scalarField& rhoM = rhoM_.internalField();
 516
 517     scalarField& dsmcRhoN = dsmcRhoN_.internalField();
 518
 519     scalarField& linearKE = linearKE_.internalField();
 520
 521     scalarField& internalE = internalE_.internalField();
 522
 523     scalarField& iDof = iDof_.internalField();
 524
 525     vectorField& momentum = momentum_.internalField();
 526
 527     forAllConstIter(typename DsmcCloud<ParcelType>, *this, iter)
 528     {
 529         const ParcelType& p = iter();
 530         const label cellI = p.cell();
 531
 532         rhoN[cellI]++;
 533
 534         rhoM[cellI] += constProps(p.typeId()).mass();
 535
 536         dsmcRhoN[cellI]++;
 537
 538         linearKE[cellI] += 0.5*constProps(p.typeId()).mass()*(p.U() & p.U());
 539
 540         internalE[cellI] += p.Ei();
 541
 542         iDof[cellI] += constProps(p.typeId()).internalDegreesOfFreedom();
 543
 544         momentum[cellI] += constProps(p.typeId()).mass()*p.U();
 545     }
 546
 547     rhoN *= nParticle_/mesh().cellVolumes();
 548     rhoN_.correctBoundaryConditions();
 549
 550     rhoM *= nParticle_/mesh().cellVolumes();
 551     rhoM_.correctBoundaryConditions();
 552
 553     linearKE *= nParticle_/mesh().cellVolumes();
 554     linearKE_.correctBoundaryConditions();
 555
 556     internalE *= nParticle_/mesh().cellVolumes();
 557     internalE_.correctBoundaryConditions();
 558
 559     iDof *= nParticle_/mesh().cellVolumes();
 560     iDof_.correctBoundaryConditions();
 561
 562     momentum *= nParticle_/mesh().cellVolumes();
 563     momentum_.correctBoundaryConditions();
 564 }
 565
 566
 567 // * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * * //
 568
 569 template<class ParcelType>
 570 void Foam::DsmcCloud<ParcelType>::addNewParcel
 571 (
 572     const vector& position,
 573     const vector& U,
 574     const scalar Ei,
 575     const label cellId,
 576     const label typeId
 577 )
 578 {
 579     ParcelType* pPtr = new ParcelType
 580     (
 581         *this,
 582         position,
 583         U,
 584         Ei,
 585         cellId,
 586         typeId
 587     );
 588
 589     addParticle(pPtr);
 590 }
 591
 592
 593 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 594
 595 template<class ParcelType>
 596 Foam::DsmcCloud<ParcelType>::DsmcCloud
 597 (
 598     const word& cloudName,
 599     const fvMesh& mesh,
 600     bool readFields
 601 )
 602 :
 603     Cloud<ParcelType>(mesh, cloudName, false),
 604     DsmcBaseCloud(),
 605     cloudName_(cloudName),
 606     mesh_(mesh),
 607     particleProperties_
 608     (
 609         IOobject
 610         (
 611             cloudName + "Properties",
 612             mesh_.time().constant(),
 613             mesh_,
 614             IOobject::MUST_READ,
 615             IOobject::NO_WRITE
 616         )
 617     ),
 618     typeIdList_(particleProperties_.lookup("typeIdList")),
 619     nParticle_(readScalar(particleProperties_.lookup("nEquivalentParticles"))),
 620     cellOccupancy_(mesh_.nCells()),
 621     sigmaTcRMax_
 622     (
 623         IOobject
 624         (
 625             this->name() + "SigmaTcRMax",
 626             mesh_.time().timeName(),
 627             mesh_,
 628             IOobject::MUST_READ,
 629             IOobject::AUTO_WRITE
 630         ),
 631         mesh_
 632     ),
 633     collisionSelectionRemainder_(mesh_.nCells(), 0),
 634     q_
 635     (
 636         IOobject
 637         (
 638             "q",
 639             mesh_.time().timeName(),
 640             mesh_,
 641             IOobject::MUST_READ,
 642             IOobject::AUTO_WRITE
 643         ),
 644         mesh_
 645     ),
 646     fD_
 647     (
 648         IOobject
 649         (
 650             "fD",
 651             mesh_.time().timeName(),
 652             mesh_,
 653             IOobject::MUST_READ,
 654             IOobject::AUTO_WRITE
 655         ),
 656         mesh_
 657     ),
 658     rhoN_
 659     (
 660         IOobject
 661         (
 662             "rhoN",
 663             mesh_.time().timeName(),
 664             mesh_,
 665             IOobject::MUST_READ,
 666             IOobject::AUTO_WRITE
 667         ),
 668         mesh_
 669     ),
 670     rhoM_
 671     (
 672         IOobject
 673         (
 674             "rhoM",
 675             mesh_.time().timeName(),
 676             mesh_,
 677             IOobject::MUST_READ,
 678             IOobject::AUTO_WRITE
 679         ),
 680         mesh_
 681     ),
 682     dsmcRhoN_
 683     (
 684         IOobject
 685         (
 686             "dsmcRhoN",
 687             mesh_.time().timeName(),
 688             mesh_,
 689             IOobject::MUST_READ,
 690             IOobject::AUTO_WRITE
 691         ),
 692         mesh_
 693     ),
 694     linearKE_
 695     (
 696         IOobject
 697         (
 698             "linearKE",
 699             mesh_.time().timeName(),
 700             mesh_,
 701             IOobject::MUST_READ,
 702             IOobject::AUTO_WRITE
 703         ),
 704         mesh_
 705     ),
 706     internalE_
 707     (
 708         IOobject
 709         (
 710             "internalE",
 711             mesh_.time().timeName(),
 712             mesh_,
 713             IOobject::MUST_READ,
 714             IOobject::AUTO_WRITE
 715         ),
 716         mesh_
 717     ),
 718     iDof_
 719     (
 720         IOobject
 721         (
 722             "iDof",
 723             mesh_.time().timeName(),
 724             mesh_,
 725             IOobject::MUST_READ,
 726             IOobject::AUTO_WRITE
 727         ),
 728         mesh_
 729     ),
 730     momentum_
 731     (
 732         IOobject
 733         (
 734             "momentum",
 735             mesh_.time().timeName(),
 736             mesh_,
 737             IOobject::MUST_READ,
 738             IOobject::AUTO_WRITE
 739         ),
 740         mesh_
 741     ),
 742     constProps_(),
 743     rndGen_(label(149382906) + 7183*Pstream::myProcNo()),
 744     boundaryT_
 745     (
 746         volScalarField
 747         (
 748             IOobject
 749             (
 750                 "boundaryT",
 751                 mesh_.time().timeName(),
 752                 mesh_,
 753                 IOobject::MUST_READ,
 754                 IOobject::AUTO_WRITE
 755             ),
 756             mesh_
 757         )
 758     ),
 759     boundaryU_
 760     (
 761         volVectorField
 762         (
 763             IOobject
 764             (
 765                 "boundaryU",
 766                 mesh_.time().timeName(),
 767                 mesh_,
 768                 IOobject::MUST_READ,
 769                 IOobject::AUTO_WRITE
 770             ),
 771             mesh_
 772         )
 773     ),
 774     binaryCollisionModel_
 775     (
 776         BinaryCollisionModel<DsmcCloud<ParcelType> >::New
 777         (
 778             particleProperties_,
 779             *this
 780         )
 781     ),
 782     wallInteractionModel_
 783     (
 784         WallInteractionModel<DsmcCloud<ParcelType> >::New
 785         (
 786             particleProperties_,
 787             *this
 788         )
 789     ),
 790     inflowBoundaryModel_
 791     (
 792         InflowBoundaryModel<DsmcCloud<ParcelType> >::New
 793         (
 794             particleProperties_,
 795             *this
 796         )
 797     )
 798 {
 799     buildConstProps();
 800
 801     buildCellOccupancy();
 802
 803     // Initialise the collision selection remainder to a random value between 0
 804     // and 1.
 805     forAll(collisionSelectionRemainder_, i)
 806     {
 807         collisionSelectionRemainder_[i] = rndGen_.scalar01();
 808     }
 809
 810     if (readFields)
 811     {
 812         ParcelType::readFields(*this);
 813     }
 814 }
 815
 816
 817 template<class ParcelType>
 818 Foam::DsmcCloud<ParcelType>::DsmcCloud
 819 (
 820     const word& cloudName,
 821     const fvMesh& mesh,
 822     const IOdictionary& dsmcInitialiseDict
 823 )
 824     :
 825     Cloud<ParcelType>(mesh, cloudName, false),
 826     DsmcBaseCloud(),
 827     cloudName_(cloudName),
 828     mesh_(mesh),
 829     particleProperties_
 830     (
 831         IOobject
 832         (
 833             cloudName + "Properties",
 834             mesh_.time().constant(),
 835             mesh_,
 836             IOobject::MUST_READ,
 837             IOobject::NO_WRITE
 838         )
 839     ),
 840     typeIdList_(particleProperties_.lookup("typeIdList")),
 841     nParticle_(readScalar(particleProperties_.lookup("nEquivalentParticles"))),
 842     cellOccupancy_(),
 843     sigmaTcRMax_
 844     (
 845         IOobject
 846         (
 847             this->name() + "SigmaTcRMax",
 848             mesh_.time().timeName(),
 849             mesh_,
 850             IOobject::NO_READ,
 851             IOobject::AUTO_WRITE
 852         ),
 853         mesh_,
 854         dimensionedScalar("zero",  dimensionSet(0, 3, -1, 0, 0), 0.0)
 855     ),
 856     collisionSelectionRemainder_(),
 857     q_
 858     (
 859         IOobject
 860         (
 861             this->name() + "q_",
 862             mesh_.time().timeName(),
 863             mesh_,
 864             IOobject::NO_READ,
 865             IOobject::NO_WRITE
 866         ),
 867         mesh_,
 868         dimensionedScalar("zero",  dimensionSet(1, 0, -3, 0, 0), 0.0)
 869     ),
 870     fD_
 871     (
 872         IOobject
 873         (
 874             this->name() + "fD_",
 875             mesh_.time().timeName(),
 876             mesh_,
 877             IOobject::NO_READ,
 878             IOobject::NO_WRITE
 879         ),
 880         mesh_,
 881         dimensionedVector
 882         (
 883             "zero",
 884             dimensionSet(1, -1, -2, 0, 0),
 885             vector::zero
 886         )
 887     ),
 888     rhoN_
 889     (
 890         IOobject
 891         (
 892             this->name() + "rhoN_",
 893             mesh_.time().timeName(),
 894             mesh_,
 895             IOobject::NO_READ,
 896             IOobject::NO_WRITE
 897         ),
 898         mesh_,
 899         dimensionedScalar("zero",  dimensionSet(0, -3, 0, 0, 0), VSMALL)
 900     ),
 901     rhoM_
 902     (
 903         IOobject
 904         (
 905             this->name() + "rhoM_",
 906             mesh_.time().timeName(),
 907             mesh_,
 908             IOobject::NO_READ,
 909             IOobject::NO_WRITE
 910         ),
 911         mesh_,
 912         dimensionedScalar("zero",  dimensionSet(1, -3, 0, 0, 0), VSMALL)
 913     ),
 914     dsmcRhoN_
 915     (
 916         IOobject
 917         (
 918             this->name() + "dsmcRhoN_",
 919             mesh_.time().timeName(),
 920             mesh_,
 921             IOobject::NO_READ,
 922             IOobject::NO_WRITE
 923         ),
 924         mesh_,
 925         dimensionedScalar("zero",  dimensionSet(0, -3, 0, 0, 0), 0.0)
 926     ),
 927     linearKE_
 928     (
 929         IOobject
 930         (
 931             this->name() + "linearKE_",
 932             mesh_.time().timeName(),
 933             mesh_,
 934             IOobject::NO_READ,
 935             IOobject::NO_WRITE
 936         ),
 937         mesh_,
 938         dimensionedScalar("zero",  dimensionSet(1, -1, -2, 0, 0), 0.0)
 939     ),
 940     internalE_
 941     (
 942         IOobject
 943         (
 944             this->name() + "internalE_",
 945             mesh_.time().timeName(),
 946             mesh_,
 947             IOobject::NO_READ,
 948             IOobject::NO_WRITE
 949         ),
 950         mesh_,
 951         dimensionedScalar("zero",  dimensionSet(1, -1, -2, 0, 0), 0.0)
 952     ),
 953     iDof_
 954     (
 955         IOobject
 956         (
 957             this->name() + "iDof_",
 958             mesh_.time().timeName(),
 959             mesh_,
 960             IOobject::NO_READ,
 961             IOobject::NO_WRITE
 962         ),
 963         mesh_,
 964         dimensionedScalar("zero",  dimensionSet(0, -3, 0, 0, 0), VSMALL)
 965     ),
 966     momentum_
 967     (
 968         IOobject
 969         (
 970             this->name() + "momentum_",
 971             mesh_.time().timeName(),
 972             mesh_,
 973             IOobject::NO_READ,
 974             IOobject::NO_WRITE
 975         ),
 976         mesh_,
 977         dimensionedVector
 978         (
 979             "zero",
 980             dimensionSet(1, -2, -1, 0, 0),
 981             vector::zero
 982         )
 983     ),
 984     constProps_(),
 985     rndGen_(label(971501) + 1526*Pstream::myProcNo()),
 986     boundaryT_
 987     (
 988         volScalarField
 989         (
 990             IOobject
 991             (
 992                 "boundaryT",
 993                 mesh_.time().timeName(),
 994                 mesh_,
 995                 IOobject::NO_READ,
 996                 IOobject::NO_WRITE
 997             ),
 998             mesh_,
 999             dimensionedScalar("zero",  dimensionSet(0, 0, 0, 1, 0), 0.0)
1000         )
1001     ),
1002     boundaryU_
1003     (
1004         volVectorField
1005         (
1006             IOobject
1007             (
1008                 "boundaryU",
1009                 mesh_.time().timeName(),
1010                 mesh_,
1011                 IOobject::NO_READ,
1012                 IOobject::NO_WRITE
1013             ),
1014             mesh_,
1015             dimensionedVector
1016             (
1017                 "zero",
1018                 dimensionSet(0, 1, -1, 0, 0),
1019                 vector::zero
1020             )
1021         )
1022     ),
1023     binaryCollisionModel_(),
1024     wallInteractionModel_(),
1025     inflowBoundaryModel_()
1026 {
1027     clear();
1028
1029     buildConstProps();
1030
1031     initialise(dsmcInitialiseDict);
1032 }
1033
1034
1035 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
1036
1037 template<class ParcelType>
1038 Foam::DsmcCloud<ParcelType>::~DsmcCloud()
1039 {}
1040
1041
1042 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
1043
1044 template<class ParcelType>
1045 void Foam::DsmcCloud<ParcelType>::evolve()
1046 {
1047     typename ParcelType::trackData td(*this);
1048
1049     // Reset the data collection fields
1050     resetFields();
1051
1052     if (debug)
1053     {
1054         this->dumpParticlePositions();
1055     }
1056
1057     // Insert new particles from the inflow boundary
1058     this->inflowBoundary().inflow();
1059
1060     // Move the particles ballistically with their current velocities
1061     Cloud<ParcelType>::move(td);
1062
1063     // Calculate new velocities via stochastic collisions
1064     collisions();
1065
1066     // Calculate the volume field data
1067     calculateFields();
1068 }
1069
1070
1071 template<class ParcelType>
1072 void Foam::DsmcCloud<ParcelType>::info() const
1073 {
1074     label nDsmcParticles = this->size();
1075     reduce(nDsmcParticles, sumOp<label>());
1076
1077     scalar nMol = nDsmcParticles*nParticle_;
1078
1079     vector linearMomentum = linearMomentumOfSystem();
1080     reduce(linearMomentum, sumOp<vector>());
1081
1082     scalar linearKineticEnergy = linearKineticEnergyOfSystem();
1083     reduce(linearKineticEnergy, sumOp<scalar>());
1084
1085     scalar internalEnergy = internalEnergyOfSystem();
1086     reduce(internalEnergy, sumOp<scalar>());
1087
1088     Info<< "Cloud name: " << this->name() << nl
1089         << "    Number of dsmc particles        = "
1090         << nDsmcParticles
1091         << endl;
1092
1093     if (nDsmcParticles)
1094     {
1095         Info<< "    Number of molecules             = "
1096             << nMol << nl
1097             << "    Mass in system                  = "
1098             << returnReduce(massInSystem(), sumOp<scalar>()) << nl
1099             << "    Average linear momentum         = "
1100             << linearMomentum/nMol << nl
1101             << "   |Average linear momentum|        = "
1102             << mag(linearMomentum)/nMol << nl
1103             << "    Average linear kinetic energy   = "
1104             << linearKineticEnergy/nMol << nl
1105             << "    Average internal energy         = "
1106             << internalEnergy/nMol << nl
1107             << "    Average total energy            = "
1108             << (internalEnergy + linearKineticEnergy)/nMol
1109             << endl;
1110     }
1111 }
1112
1113
1114 template<class ParcelType>
1115 Foam::vector Foam::DsmcCloud<ParcelType>::equipartitionLinearVelocity
1116 (
1117     scalar temperature,
1118     scalar mass
1119 )
1120 {
1121     return
1122         sqrt(kb*temperature/mass)
1123        *vector
1124         (
1125             rndGen_.GaussNormal(),
1126             rndGen_.GaussNormal(),
1127             rndGen_.GaussNormal()
1128         );
1129 }
1130
1131
1132 template<class ParcelType>
1133 Foam::scalar Foam::DsmcCloud<ParcelType>::equipartitionInternalEnergy
1134 (
1135     scalar temperature,
1136     scalar iDof
1137 )
1138 {
1139     scalar Ei = 0.0;
1140
1141     if (iDof < SMALL)
1142     {
1143         return Ei;
1144     }
1145     else if (iDof < 2.0 + SMALL && iDof > 2.0 - SMALL)
1146     {
1147         // Special case for iDof = 2, i.e. diatomics;
1148         Ei = -log(rndGen_.scalar01())*kb*temperature;
1149     }
1150     else
1151     {
1152         scalar a = 0.5*iDof - 1;
1153
1154         scalar energyRatio;
1155
1156         scalar P = -1;
1157
1158         do
1159         {
1160             energyRatio = 10*rndGen_.scalar01();
1161
1162             P = pow((energyRatio/a), a)*exp(a - energyRatio);
1163
1164         } while (P < rndGen_.scalar01());
1165
1166         Ei = energyRatio*kb*temperature;
1167     }
1168
1169     return Ei;
1170 }
1171
1172
1173 template<class ParcelType>
1174 void Foam::DsmcCloud<ParcelType>::dumpParticlePositions() const
1175 {
1176     OFstream pObj
1177     (
1178         this->db().time().path()/"parcelPositions_"
1179       + this->name() + "_"
1180       + this->db().time().timeName() + ".obj"
1181     );
1182
1183     forAllConstIter(typename DsmcCloud<ParcelType>, *this, iter)
1184     {
1185         const ParcelType& p = iter();
1186
1187         pObj<< "v " << p.position().x()
1188             << " "  << p.position().y()
1189             << " "  << p.position().z()
1190             << nl;
1191     }
1192
1193     pObj.flush();
1194 }
1195
1196
1197 // ************************************************************************* //