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 = cachedDeltaT();
 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>::resetSurfaceDataFields()
 479 {
 480     volScalarField::GeometricBoundaryField& qBF(q_.boundaryField());
 481
 482     forAll(qBF, p)
 483     {
 484         qBF[p] = 0.0;
 485     }
 486
 487     volVectorField::GeometricBoundaryField& fDBF(fD_.boundaryField());
 488
 489     forAll(fDBF, p)
 490     {
 491         fDBF[p] = vector::zero;
 492     }
 493 }
 494
 495
 496 // * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * * //
 497
 498 template<class ParcelType>
 499 void Foam::DsmcCloud<ParcelType>::addNewParcel
 500 (
 501     const vector& position,
 502     const vector& U,
 503     const scalar Ei,
 504     const label cellId,
 505     const label typeId
 506 )
 507 {
 508     ParcelType* pPtr = new ParcelType
 509     (
 510         *this,
 511         position,
 512         U,
 513         Ei,
 514         cellId,
 515         typeId
 516     );
 517
 518     addParticle(pPtr);
 519 }
 520
 521
 522 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 523
 524 template<class ParcelType>
 525 Foam::DsmcCloud<ParcelType>::DsmcCloud
 526 (
 527     const word& cloudName,
 528     const volScalarField& T,
 529     const volVectorField& U
 530 )
 531 :
 532     Cloud<ParcelType>(T.mesh(), cloudName, false),
 533     DsmcBaseCloud(),
 534     cloudName_(cloudName),
 535     mesh_(T.mesh()),
 536     particleProperties_
 537     (
 538         IOobject
 539         (
 540             cloudName + "Properties",
 541             mesh_.time().constant(),
 542             mesh_,
 543             IOobject::MUST_READ,
 544             IOobject::NO_WRITE
 545         )
 546     ),
 547     typeIdList_(particleProperties_.lookup("typeIdList")),
 548     nParticle_(readScalar(particleProperties_.lookup("nEquivalentParticles"))),
 549     cellOccupancy_(mesh_.nCells()),
 550     sigmaTcRMax_
 551     (
 552         IOobject
 553         (
 554             this->name() + "SigmaTcRMax",
 555             mesh_.time().timeName(),
 556             mesh_,
 557             IOobject::MUST_READ,
 558             IOobject::AUTO_WRITE
 559         ),
 560         mesh_
 561     ),
 562     collisionSelectionRemainder_(mesh_.nCells(), 0),
 563     q_
 564     (
 565         IOobject
 566         (
 567             this->name() + "q_",
 568             mesh_.time().timeName(),
 569             mesh_,
 570             IOobject::NO_READ,
 571             IOobject::NO_WRITE
 572         ),
 573         mesh_,
 574         dimensionedScalar("zero",  dimensionSet(1, 0, -3, 0, 0), 0.0)
 575     ),
 576     fD_
 577     (
 578         IOobject
 579         (
 580             this->name() + "fD_",
 581             mesh_.time().timeName(),
 582             mesh_,
 583             IOobject::NO_READ,
 584             IOobject::NO_WRITE
 585         ),
 586         mesh_,
 587         dimensionedVector
 588         (
 589             "zero",
 590             dimensionSet(1, -1, -2, 0, 0),
 591             vector::zero
 592         )
 593     ),
 594     constProps_(),
 595     rndGen_(label(149382906) + 7183*Pstream::myProcNo()),
 596     T_(T),
 597     U_(U),
 598     binaryCollisionModel_
 599     (
 600         BinaryCollisionModel<DsmcCloud<ParcelType> >::New
 601         (
 602             particleProperties_,
 603             *this
 604         )
 605     ),
 606     wallInteractionModel_
 607     (
 608         WallInteractionModel<DsmcCloud<ParcelType> >::New
 609         (
 610             particleProperties_,
 611             *this
 612         )
 613     ),
 614     inflowBoundaryModel_
 615     (
 616         InflowBoundaryModel<DsmcCloud<ParcelType> >::New
 617         (
 618             particleProperties_,
 619             *this
 620         )
 621     )
 622 {
 623     buildConstProps();
 624
 625     buildCellOccupancy();
 626
 627     // Initialise the collision selection remainder to a random value between 0
 628     // and 1.
 629     forAll(collisionSelectionRemainder_, i)
 630     {
 631         collisionSelectionRemainder_[i] = rndGen_.scalar01();
 632     }
 633 }
 634
 635
 636 template<class ParcelType>
 637 Foam::DsmcCloud<ParcelType>::DsmcCloud
 638 (
 639     const word& cloudName,
 640     const fvMesh& mesh
 641 )
 642     :
 643     Cloud<ParcelType>(mesh, cloudName, false),
 644     DsmcBaseCloud(),
 645     cloudName_(cloudName),
 646     mesh_(mesh),
 647     particleProperties_
 648     (
 649         IOobject
 650         (
 651             cloudName + "Properties",
 652             mesh_.time().constant(),
 653             mesh_,
 654             IOobject::MUST_READ,
 655             IOobject::NO_WRITE
 656         )
 657     ),
 658     typeIdList_(particleProperties_.lookup("typeIdList")),
 659     nParticle_(readScalar(particleProperties_.lookup("nEquivalentParticles"))),
 660     cellOccupancy_(),
 661     sigmaTcRMax_
 662     (
 663         IOobject
 664         (
 665             this->name() + "SigmaTcRMax",
 666             mesh_.time().timeName(),
 667             mesh_,
 668             IOobject::NO_READ,
 669             IOobject::AUTO_WRITE
 670         ),
 671         mesh_,
 672         dimensionedScalar("zero",  dimensionSet(0, 3, -1, 0, 0), 0.0)
 673     ),
 674     collisionSelectionRemainder_(),
 675     q_
 676     (
 677         IOobject
 678         (
 679             this->name() + "q_",
 680             mesh_.time().timeName(),
 681             mesh_,
 682             IOobject::NO_READ,
 683             IOobject::NO_WRITE
 684         ),
 685         mesh_,
 686         dimensionedScalar("zero",  dimensionSet(1, 0, -3, 0, 0), 0.0)
 687     ),
 688     fD_
 689     (
 690         IOobject
 691         (
 692             this->name() + "fD_",
 693             mesh_.time().timeName(),
 694             mesh_,
 695             IOobject::NO_READ,
 696             IOobject::NO_WRITE
 697         ),
 698         mesh_,
 699         dimensionedVector
 700         (
 701             "zero",
 702             dimensionSet(1, -1, -2, 0, 0),
 703             vector::zero
 704         )
 705     ),
 706     constProps_(),
 707     rndGen_(label(971501) + 1526*Pstream::myProcNo()),
 708     T_
 709     (
 710         volScalarField
 711         (
 712             IOobject
 713             (
 714                 "T",
 715                 mesh_.time().timeName(),
 716                 mesh_,
 717                 IOobject::NO_READ,
 718                 IOobject::NO_WRITE
 719             ),
 720             mesh_,
 721             dimensionedScalar("zero",  dimensionSet(0, 0, 0, 1, 0), 0.0)
 722         )
 723     ),
 724     U_
 725     (
 726         volVectorField
 727         (
 728             IOobject
 729             (
 730                 "U",
 731                 mesh_.time().timeName(),
 732                 mesh_,
 733                 IOobject::NO_READ,
 734                 IOobject::NO_WRITE
 735             ),
 736             mesh_,
 737             dimensionedVector
 738             (
 739                 "zero",
 740                 dimensionSet(0, 1, -1, 0, 0),
 741                 vector::zero
 742             )
 743         )
 744     ),
 745     binaryCollisionModel_(),
 746     wallInteractionModel_(),
 747     inflowBoundaryModel_()
 748 {
 749     clear();
 750
 751     buildConstProps();
 752
 753     IOdictionary dsmcInitialiseDict
 754     (
 755         IOobject
 756         (
 757             "dsmcInitialiseDict",
 758             mesh_.time().system(),
 759             mesh_,
 760             IOobject::MUST_READ,
 761             IOobject::NO_WRITE
 762         )
 763     );
 764
 765     initialise(dsmcInitialiseDict);
 766 }
 767
 768
 769 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 770
 771 template<class ParcelType>
 772 Foam::DsmcCloud<ParcelType>::~DsmcCloud()
 773 {}
 774
 775
 776 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 777
 778 template<class ParcelType>
 779 void Foam::DsmcCloud<ParcelType>::evolve()
 780 {
 781     // cache the value of deltaT for this timestep
 782     storeDeltaT();
 783
 784     typename ParcelType::trackData td(*this);
 785
 786     // Reset the surface data collection fields
 787     resetSurfaceDataFields();
 788
 789     if (debug)
 790     {
 791         this->dumpParticlePositions();
 792     }
 793
 794     // Insert new particles from the inflow boundary
 795     this->inflowBoundary().inflow();
 796
 797     // Move the particles ballistically with their current velocities
 798     Cloud<ParcelType>::move(td);
 799
 800     // Calculate new velocities via stochastic collisions
 801     collisions();
 802 }
 803
 804
 805 template<class ParcelType>
 806 void Foam::DsmcCloud<ParcelType>::info() const
 807 {
 808     label nDsmcParticles = this->size();
 809     reduce(nDsmcParticles, sumOp<label>());
 810
 811     scalar nMol = nDsmcParticles*nParticle_;
 812
 813     vector linearMomentum = linearMomentumOfSystem();
 814     reduce(linearMomentum, sumOp<vector>());
 815
 816     scalar linearKineticEnergy = linearKineticEnergyOfSystem();
 817     reduce(linearKineticEnergy, sumOp<scalar>());
 818
 819     scalar internalEnergy = internalEnergyOfSystem();
 820     reduce(internalEnergy, sumOp<scalar>());
 821
 822     Info<< "Cloud name: " << this->name() << nl
 823         << "    Number of dsmc particles        = "
 824         << nDsmcParticles
 825         << endl;
 826
 827     if (nDsmcParticles)
 828     {
 829         Info<< "    Number of molecules             = "
 830             << nMol << nl
 831             << "    Mass in system                  = "
 832             << returnReduce(massInSystem(), sumOp<scalar>()) << nl
 833             << "    Average linear momentum         = "
 834             << linearMomentum/nMol << nl
 835             << "   |Average linear momentum|        = "
 836             << mag(linearMomentum)/nMol << nl
 837             << "    Average linear kinetic energy   = "
 838             << linearKineticEnergy/nMol << nl
 839             << "    Average internal energy         = "
 840             << internalEnergy/nMol << nl
 841             << "    Average total energy            = "
 842             << (internalEnergy + linearKineticEnergy)/nMol
 843             << endl;
 844     }
 845 }
 846
 847
 848 template<class ParcelType>
 849 Foam::vector Foam::DsmcCloud<ParcelType>::equipartitionLinearVelocity
 850 (
 851     scalar temperature,
 852     scalar mass
 853 )
 854 {
 855     return
 856         sqrt(kb*temperature/mass)
 857        *vector
 858         (
 859             rndGen_.GaussNormal(),
 860             rndGen_.GaussNormal(),
 861             rndGen_.GaussNormal()
 862         );
 863 }
 864
 865
 866 template<class ParcelType>
 867 Foam::scalar Foam::DsmcCloud<ParcelType>::equipartitionInternalEnergy
 868 (
 869     scalar temperature,
 870     scalar iDof
 871 )
 872 {
 873     scalar Ei = 0.0;
 874
 875     if (iDof < SMALL)
 876     {
 877         return Ei;
 878     }
 879     else if (iDof < 2.0 + SMALL && iDof > 2.0 - SMALL)
 880     {
 881         // Special case for iDof = 2, i.e. diatomics;
 882         Ei = -log(rndGen_.scalar01())*kb*temperature;
 883     }
 884     else
 885     {
 886         scalar a = 0.5*iDof - 1;
 887
 888         scalar energyRatio;
 889
 890         scalar P = -1;
 891
 892         do
 893         {
 894             energyRatio = 10*rndGen_.scalar01();
 895
 896             P = pow((energyRatio/a), a)*exp(a - energyRatio);
 897
 898         } while (P < rndGen_.scalar01());
 899
 900         Ei = energyRatio*kb*temperature;
 901     }
 902
 903     return Ei;
 904 }
 905
 906
 907 template<class ParcelType>
 908 void Foam::DsmcCloud<ParcelType>::dumpParticlePositions() const
 909 {
 910     OFstream pObj
 911     (
 912         this->db().time().path()/"parcelPositions_"
 913       + this->name() + "_"
 914       + this->db().time().timeName() + ".obj"
 915     );
 916
 917     forAllConstIter(typename DsmcCloud<ParcelType>, *this, iter)
 918     {
 919         const ParcelType& p = iter();
 920
 921         pObj<< "v " << p.position().x()
 922             << " "  << p.position().y()
 923             << " "  << p.position().z()
 924             << nl;
 925     }
 926
 927     pObj.flush();
 928 }
 929
 930
 931 // ************************************************************************* //