DSMCCloud.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2020 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 "DSMCCloud.H"
27 #include "NoBinaryCollision.H"
28 #include "WallInteractionModel.H"
29 #include "InflowBoundaryModel.H"
30 #include "constants.H"
33 
34 using namespace Foam::constant;
35 
36 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
37 
38 template<class ParcelType>
40 {
41  Info<< nl << "Constructing constant properties for" << endl;
42  constProps_.setSize(typeIdList_.size());
43 
44  dictionary moleculeProperties
45  (
46  particleProperties_.subDict("moleculeProperties")
47  );
48 
49  forAll(typeIdList_, i)
50  {
51  const word& id(typeIdList_[i]);
52 
53  Info<< " " << id << endl;
54 
55  const dictionary& molDict(moleculeProperties.subDict(id));
56 
57  constProps_[i] =
58  typename ParcelType::constantProperties(molDict);
59  }
60 }
61 
62 
63 template<class ParcelType>
65 {
66  forAll(cellOccupancy_, cO)
67  {
68  cellOccupancy_[cO].clear();
69  }
70 
71  forAllIter(typename DSMCCloud<ParcelType>, *this, iter)
72  {
73  cellOccupancy_[iter().cell()].append(&iter());
74  }
75 }
76 
77 
78 template<class ParcelType>
80 (
81  const IOdictionary& dsmcInitialiseDict
82 )
83 {
84  Info<< nl << "Initialising particles" << endl;
85 
86  const scalar temperature
87  (
88  dsmcInitialiseDict.template lookup<scalar>("temperature")
89  );
90 
91  const vector velocity(dsmcInitialiseDict.lookup("velocity"));
92 
93  const dictionary& numberDensitiesDict
94  (
95  dsmcInitialiseDict.subDict("numberDensities")
96  );
97 
98  List<word> molecules(numberDensitiesDict.toc());
99 
100  Field<scalar> numberDensities(molecules.size());
101 
102  forAll(molecules, i)
103  {
104  numberDensities[i] =
105  numberDensitiesDict.lookup<scalar>(molecules[i]);
106  }
107 
108  numberDensities /= nParticle_;
109 
110  forAll(mesh_.cells(), celli)
111  {
112  List<tetIndices> cellTets = polyMeshTetDecomposition::cellTetIndices
113  (
114  mesh_,
115  celli
116  );
117 
118  forAll(cellTets, tetI)
119  {
120  const tetIndices& cellTetIs = cellTets[tetI];
121  tetPointRef tet = cellTetIs.tet(mesh_);
122  scalar tetVolume = tet.mag();
123 
124  forAll(molecules, i)
125  {
126  const word& moleculeName(molecules[i]);
127 
128  label typeId(findIndex(typeIdList_, moleculeName));
129 
130  if (typeId == -1)
131  {
133  << "typeId " << moleculeName << "not defined." << nl
134  << abort(FatalError);
135  }
136 
137  const typename ParcelType::constantProperties& cP =
138  constProps(typeId);
139 
140  scalar numberDensity = numberDensities[i];
141 
142  // Calculate the number of particles required
143  scalar particlesRequired = numberDensity*tetVolume;
144 
145  // Only integer numbers of particles can be inserted
146  label nParticlesToInsert = label(particlesRequired);
147 
148  // Add another particle with a probability proportional to the
149  // remainder of taking the integer part of particlesRequired
150  if
151  (
152  (particlesRequired - nParticlesToInsert)
153  > rndGen_.scalar01()
154  )
155  {
156  nParticlesToInsert++;
157  }
158 
159  for (label pI = 0; pI < nParticlesToInsert; pI++)
160  {
161  point p = tet.randomPoint(rndGen_);
162 
163  vector U = equipartitionLinearVelocity
164  (
165  temperature,
166  cP.mass()
167  );
168 
169  scalar Ei = equipartitionInternalEnergy
170  (
171  temperature,
172  cP.internalDegreesOfFreedom()
173  );
174 
175  U += velocity;
176 
177  addNewParcel(p, celli, U, Ei, typeId);
178  }
179  }
180  }
181  }
182 
183  // Initialise the sigmaTcRMax_ field to the product of the cross section of
184  // the most abundant species and the most probable thermal speed (Bird,
185  // p222-223)
186 
187  label mostAbundantType(findMax(numberDensities));
188 
189  const typename ParcelType::constantProperties& cP = constProps
190  (
191  mostAbundantType
192  );
193 
194  sigmaTcRMax_.primitiveFieldRef() = cP.sigmaT()*maxwellianMostProbableSpeed
195  (
196  temperature,
197  cP.mass()
198  );
199 
200  sigmaTcRMax_.correctBoundaryConditions();
201 }
202 
203 
204 template<class ParcelType>
206 {
207  if (isType<NoBinaryCollision<DSMCCloud<ParcelType>>>(binaryCollision()))
208  {
209  return;
210  }
211 
212  // Temporary storage for subCells
213  List<DynamicList<label>> subCells(8);
214 
215  scalar deltaT = mesh().time().deltaTValue();
216 
217  label collisionCandidates = 0;
218 
219  label collisions = 0;
220 
221  forAll(cellOccupancy_, celli)
222  {
223  const DynamicList<ParcelType*>& cellParcels(cellOccupancy_[celli]);
224 
225  label nC(cellParcels.size());
226 
227  if (nC > 1)
228  {
229  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
230  // Assign particles to one of 8 Cartesian subCells
231 
232  // Clear temporary lists
233  forAll(subCells, i)
234  {
235  subCells[i].clear();
236  }
237 
238  // Inverse addressing specifying which subCell a parcel is in
239  List<label> whichSubCell(cellParcels.size());
240 
241  const point& cC = mesh_.cellCentres()[celli];
242 
243  forAll(cellParcels, i)
244  {
245  const ParcelType& p = *cellParcels[i];
246  vector relPos = p.position() - cC;
247 
248  label subCell =
249  pos0(relPos.x()) + 2*pos0(relPos.y()) + 4*pos0(relPos.z());
250 
251  subCells[subCell].append(i);
252  whichSubCell[i] = subCell;
253  }
254 
255  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
256 
257  scalar sigmaTcRMax = sigmaTcRMax_[celli];
258 
259  scalar selectedPairs =
260  collisionSelectionRemainder_[celli]
261  + 0.5*nC*(nC - 1)*nParticle_*sigmaTcRMax*deltaT
262  /mesh_.cellVolumes()[celli];
263 
264  label nCandidates(selectedPairs);
265  collisionSelectionRemainder_[celli] = selectedPairs - nCandidates;
266  collisionCandidates += nCandidates;
267 
268  for (label c = 0; c < nCandidates; c++)
269  {
270  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
271  // subCell candidate selection procedure
272 
273  // Select the first collision candidate
274  label candidateP = rndGen_.sampleAB<label>(0, nC);
275 
276  // Declare the second collision candidate
277  label candidateQ = -1;
278 
279  List<label> subCellPs = subCells[whichSubCell[candidateP]];
280  label nSC = subCellPs.size();
281 
282  if (nSC > 1)
283  {
284  // If there are two or more particle in a subCell, choose
285  // another from the same cell. If the same candidate is
286  // chosen, choose again.
287 
288  do
289  {
290  candidateQ = subCellPs[rndGen_.sampleAB<label>(0, nSC)];
291  } while (candidateP == candidateQ);
292  }
293  else
294  {
295  // Select a possible second collision candidate from the
296  // whole cell. If the same candidate is chosen, choose
297  // again.
298 
299  do
300  {
301  candidateQ = rndGen_.sampleAB<label>(0, nC);
302  } while (candidateP == candidateQ);
303  }
304 
305  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
306  // uniform candidate selection procedure
307 
308  // // Select the first collision candidate
309  // label candidateP = rndGen_.sampleAB<label>(0, nC);
310 
311  // // Select a possible second collision candidate
312  // label candidateQ = rndGen_.sampleAB<label>(0, nC);
313 
314  // // If the same candidate is chosen, choose again
315  // while (candidateP == candidateQ)
316  // {
317  // candidateQ = rndGen_.sampleAB<label>(0, nC);
318  // }
319 
320  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
321 
322  ParcelType& parcelP = *cellParcels[candidateP];
323  ParcelType& parcelQ = *cellParcels[candidateQ];
324 
325  scalar sigmaTcR = binaryCollision().sigmaTcR
326  (
327  parcelP,
328  parcelQ
329  );
330 
331  // Update the maximum value of sigmaTcR stored, but use the
332  // initial value in the acceptance-rejection criteria because
333  // the number of collision candidates selected was based on this
334 
335  if (sigmaTcR > sigmaTcRMax_[celli])
336  {
337  sigmaTcRMax_[celli] = sigmaTcR;
338  }
339 
340  if ((sigmaTcR/sigmaTcRMax) > rndGen_.scalar01())
341  {
342  binaryCollision().collide
343  (
344  parcelP,
345  parcelQ
346  );
347 
348  collisions++;
349  }
350  }
351  }
352  }
353 
354  reduce(collisions, sumOp<label>());
355 
356  reduce(collisionCandidates, sumOp<label>());
357 
358  sigmaTcRMax_.correctBoundaryConditions();
359 
360  if (collisionCandidates)
361  {
362  Info<< " Collisions = "
363  << collisions << nl
364  << " Acceptance rate = "
365  << scalar(collisions)/scalar(collisionCandidates) << nl
366  << endl;
367  }
368  else
369  {
370  Info<< " No collisions" << endl;
371  }
372 }
373 
374 
375 template<class ParcelType>
377 {
378  q_ = dimensionedScalar( dimensionSet(1, 0, -3, 0, 0), 0);
379 
380  fD_ = dimensionedVector
381  (
382  "zero",
383  dimensionSet(1, -1, -2, 0, 0),
384  Zero
385  );
386 
387  rhoN_ = dimensionedScalar( dimensionSet(0, -3, 0, 0, 0), vSmall);
388  rhoM_ = dimensionedScalar( dimensionSet(1, -3, 0, 0, 0), vSmall);
389  dsmcRhoN_ = dimensionedScalar( dimensionSet(0, -3, 0, 0, 0), 0);
390  linearKE_ = dimensionedScalar( dimensionSet(1, -1, -2, 0, 0), 0);
391  internalE_ = dimensionedScalar( dimensionSet(1, -1, -2, 0, 0), 0);
392  iDof_ = dimensionedScalar( dimensionSet(0, -3, 0, 0, 0), vSmall);
393 
394  momentum_ = dimensionedVector
395  (
396  "zero",
397  dimensionSet(1, -2, -1, 0, 0),
398  Zero
399  );
400 }
401 
402 
403 template<class ParcelType>
405 {
406  scalarField& rhoN = rhoN_.primitiveFieldRef();
407  scalarField& rhoM = rhoM_.primitiveFieldRef();
408  scalarField& dsmcRhoN = dsmcRhoN_.primitiveFieldRef();
409  scalarField& linearKE = linearKE_.primitiveFieldRef();
410  scalarField& internalE = internalE_.primitiveFieldRef();
411  scalarField& iDof = iDof_.primitiveFieldRef();
412  vectorField& momentum = momentum_.primitiveFieldRef();
413 
414  forAllConstIter(typename DSMCCloud<ParcelType>, *this, iter)
415  {
416  const ParcelType& p = iter();
417  const label celli = p.cell();
418 
419  rhoN[celli]++;
420  rhoM[celli] += constProps(p.typeId()).mass();
421  dsmcRhoN[celli]++;
422  linearKE[celli] += 0.5*constProps(p.typeId()).mass()*(p.U() & p.U());
423  internalE[celli] += p.Ei();
424  iDof[celli] += constProps(p.typeId()).internalDegreesOfFreedom();
425  momentum[celli] += constProps(p.typeId()).mass()*p.U();
426  }
427 
428  rhoN *= nParticle_/mesh().cellVolumes();
429  rhoN_.correctBoundaryConditions();
430 
431  rhoM *= nParticle_/mesh().cellVolumes();
432  rhoM_.correctBoundaryConditions();
433 
434  dsmcRhoN_.correctBoundaryConditions();
435 
436  linearKE *= nParticle_/mesh().cellVolumes();
437  linearKE_.correctBoundaryConditions();
438 
439  internalE *= nParticle_/mesh().cellVolumes();
440  internalE_.correctBoundaryConditions();
441 
442  iDof *= nParticle_/mesh().cellVolumes();
443  iDof_.correctBoundaryConditions();
444 
445  momentum *= nParticle_/mesh().cellVolumes();
446  momentum_.correctBoundaryConditions();
447 }
448 
449 
450 // * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * * //
451 
452 template<class ParcelType>
454 (
455  const vector& position,
456  const label celli,
457  const vector& U,
458  const scalar Ei,
459  const label typeId
460 )
461 {
462  this->addParticle(new ParcelType(mesh_, position, celli, U, Ei, typeId));
463 }
464 
465 
466 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
467 
468 template<class ParcelType>
470 (
471  const word& cloudName,
472  const fvMesh& mesh,
473  bool readFields
474 )
475 :
477  cloudName_(cloudName),
478  mesh_(mesh),
479  particleProperties_
480  (
481  IOobject
482  (
483  cloudName + "Properties",
484  mesh_.time().constant(),
485  mesh_,
488  )
489  ),
490  typeIdList_(particleProperties_.lookup("typeIdList")),
491  nParticle_
492  (
493  particleProperties_.template lookup<scalar>("nEquivalentParticles")
494  ),
495  cellOccupancy_(mesh_.nCells()),
496  sigmaTcRMax_
497  (
498  IOobject
499  (
500  this->name() + "SigmaTcRMax",
501  mesh_.time().timeName(),
502  mesh_,
505  ),
506  mesh_
507  ),
508  collisionSelectionRemainder_
509  (
510  IOobject
511  (
512  this->name() + ":collisionSelectionRemainder",
513  mesh_.time().timeName(),
514  mesh_
515  ),
516  mesh_,
518  ),
519  q_
520  (
521  IOobject
522  (
523  "q",
524  mesh_.time().timeName(),
525  mesh_,
528  ),
529  mesh_
530  ),
531  fD_
532  (
533  IOobject
534  (
535  "fD",
536  mesh_.time().timeName(),
537  mesh_,
540  ),
541  mesh_
542  ),
543  rhoN_
544  (
545  IOobject
546  (
547  "rhoN",
548  mesh_.time().timeName(),
549  mesh_,
552  ),
553  mesh_
554  ),
555  rhoM_
556  (
557  IOobject
558  (
559  "rhoM",
560  mesh_.time().timeName(),
561  mesh_,
564  ),
565  mesh_
566  ),
567  dsmcRhoN_
568  (
569  IOobject
570  (
571  "dsmcRhoN",
572  mesh_.time().timeName(),
573  mesh_,
576  ),
577  mesh_
578  ),
579  linearKE_
580  (
581  IOobject
582  (
583  "linearKE",
584  mesh_.time().timeName(),
585  mesh_,
588  ),
589  mesh_
590  ),
591  internalE_
592  (
593  IOobject
594  (
595  "internalE",
596  mesh_.time().timeName(),
597  mesh_,
600  ),
601  mesh_
602  ),
603  iDof_
604  (
605  IOobject
606  (
607  "iDof",
608  mesh_.time().timeName(),
609  mesh_,
612  ),
613  mesh_
614  ),
615  momentum_
616  (
617  IOobject
618  (
619  "momentum",
620  mesh_.time().timeName(),
621  mesh_,
624  ),
625  mesh_
626  ),
627  constProps_(),
628  rndGen_(label(149382906) + 7183*Pstream::myProcNo()),
629  boundaryT_
630  (
632  (
633  IOobject
634  (
635  "boundaryT",
636  mesh_.time().timeName(),
637  mesh_,
640  ),
641  mesh_
642  )
643  ),
644  boundaryU_
645  (
647  (
648  IOobject
649  (
650  "boundaryU",
651  mesh_.time().timeName(),
652  mesh_,
655  ),
656  mesh_
657  )
658  ),
659  binaryCollisionModel_
660  (
662  (
663  particleProperties_,
664  *this
665  )
666  ),
667  wallInteractionModel_
668  (
670  (
671  particleProperties_,
672  *this
673  )
674  ),
675  inflowBoundaryModel_
676  (
678  (
679  particleProperties_,
680  *this
681  )
682  )
683 {
684  buildConstProps();
685  buildCellOccupancy();
686 
687  // Initialise the collision selection remainder to a random value between 0
688  // and 1.
689  forAll(collisionSelectionRemainder_, i)
690  {
691  collisionSelectionRemainder_[i] = rndGen_.scalar01();
692  }
693 
694  if (readFields)
695  {
696  ParcelType::readFields(*this);
697  }
698 }
699 
700 
701 template<class ParcelType>
703 (
704  const word& cloudName,
705  const fvMesh& mesh,
706  const IOdictionary& dsmcInitialiseDict
707 )
708 :
710  cloudName_(cloudName),
711  mesh_(mesh),
712  particleProperties_
713  (
714  IOobject
715  (
716  cloudName + "Properties",
717  mesh_.time().constant(),
718  mesh_,
721  )
722  ),
723  typeIdList_(particleProperties_.lookup("typeIdList")),
724  nParticle_
725  (
726  particleProperties_.template lookup<scalar>("nEquivalentParticles")
727  ),
728  cellOccupancy_(),
729  sigmaTcRMax_
730  (
731  IOobject
732  (
733  this->name() + "SigmaTcRMax",
734  mesh_.time().timeName(),
735  mesh_,
738  ),
739  mesh_,
740  dimensionedScalar( dimensionSet(0, 3, -1, 0, 0), 0),
741  zeroGradientFvPatchScalarField::typeName
742  ),
743  collisionSelectionRemainder_
744  (
745  IOobject
746  (
747  this->name() + ":collisionSelectionRemainder",
748  mesh_.time().timeName(),
749  mesh_
750  ),
751  mesh_,
753  ),
754  q_
755  (
756  IOobject
757  (
758  this->name() + "q_",
759  mesh_.time().timeName(),
760  mesh_,
763  ),
764  mesh_,
765  dimensionedScalar( dimensionSet(1, 0, -3, 0, 0), 0)
766  ),
767  fD_
768  (
769  IOobject
770  (
771  this->name() + "fD_",
772  mesh_.time().timeName(),
773  mesh_,
776  ),
777  mesh_,
779  (
780  "zero",
781  dimensionSet(1, -1, -2, 0, 0),
782  Zero
783  )
784  ),
785  rhoN_
786  (
787  IOobject
788  (
789  this->name() + "rhoN_",
790  mesh_.time().timeName(),
791  mesh_,
794  ),
795  mesh_,
796  dimensionedScalar( dimensionSet(0, -3, 0, 0, 0), vSmall)
797  ),
798  rhoM_
799  (
800  IOobject
801  (
802  this->name() + "rhoM_",
803  mesh_.time().timeName(),
804  mesh_,
807  ),
808  mesh_,
809  dimensionedScalar( dimensionSet(1, -3, 0, 0, 0), vSmall)
810  ),
811  dsmcRhoN_
812  (
813  IOobject
814  (
815  this->name() + "dsmcRhoN_",
816  mesh_.time().timeName(),
817  mesh_,
820  ),
821  mesh_,
822  dimensionedScalar( dimensionSet(0, -3, 0, 0, 0), 0)
823  ),
824  linearKE_
825  (
826  IOobject
827  (
828  this->name() + "linearKE_",
829  mesh_.time().timeName(),
830  mesh_,
833  ),
834  mesh_,
835  dimensionedScalar( dimensionSet(1, -1, -2, 0, 0), 0)
836  ),
837  internalE_
838  (
839  IOobject
840  (
841  this->name() + "internalE_",
842  mesh_.time().timeName(),
843  mesh_,
846  ),
847  mesh_,
848  dimensionedScalar( dimensionSet(1, -1, -2, 0, 0), 0)
849  ),
850  iDof_
851  (
852  IOobject
853  (
854  this->name() + "iDof_",
855  mesh_.time().timeName(),
856  mesh_,
859  ),
860  mesh_,
861  dimensionedScalar( dimensionSet(0, -3, 0, 0, 0), vSmall)
862  ),
863  momentum_
864  (
865  IOobject
866  (
867  this->name() + "momentum_",
868  mesh_.time().timeName(),
869  mesh_,
872  ),
873  mesh_,
875  (
876  "zero",
877  dimensionSet(1, -2, -1, 0, 0),
878  Zero
879  )
880  ),
881  constProps_(),
882  rndGen_(label(971501) + 1526*Pstream::myProcNo()),
883  boundaryT_
884  (
886  (
887  IOobject
888  (
889  "boundaryT",
890  mesh_.time().timeName(),
891  mesh_,
894  ),
895  mesh_,
896  dimensionedScalar( dimensionSet(0, 0, 0, 1, 0), 0)
897  )
898  ),
899  boundaryU_
900  (
902  (
903  IOobject
904  (
905  "boundaryU",
906  mesh_.time().timeName(),
907  mesh_,
910  ),
911  mesh_,
913  (
914  "zero",
915  dimensionSet(0, 1, -1, 0, 0),
916  Zero
917  )
918  )
919  ),
920  binaryCollisionModel_(),
921  wallInteractionModel_(),
922  inflowBoundaryModel_()
923 {
924  clear();
925  buildConstProps();
926  initialise(dsmcInitialiseDict);
927 }
928 
929 
930 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
931 
932 template<class ParcelType>
934 {}
935 
936 
937 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
938 
939 template<class ParcelType>
941 {
942  typename ParcelType::trackingData td(*this);
943 
944  // Reset the data collection fields
945  resetFields();
946 
947  if (debug)
948  {
949  this->dumpParticlePositions();
950  }
951 
952  // Insert new particles from the inflow boundary
953  this->inflowBoundary().inflow();
954 
955  // Move the particles ballistically with their current velocities
956  Cloud<ParcelType>::move(*this, td, mesh_.time().deltaTValue());
957 
958  // Update cell occupancy
959  buildCellOccupancy();
960 
961  // Calculate new velocities via stochastic collisions
962  collisions();
963 
964  // Calculate the volume field data
965  calculateFields();
966 }
967 
968 
969 template<class ParcelType>
971 {
972  label nDSMCParticles = this->size();
973  reduce(nDSMCParticles, sumOp<label>());
974 
975  scalar nMol = nDSMCParticles*nParticle_;
976 
977  vector linearMomentum = linearMomentumOfSystem();
978  reduce(linearMomentum, sumOp<vector>());
979 
980  scalar linearKineticEnergy = linearKineticEnergyOfSystem();
981  reduce(linearKineticEnergy, sumOp<scalar>());
982 
983  scalar internalEnergy = internalEnergyOfSystem();
984  reduce(internalEnergy, sumOp<scalar>());
985 
986  Info<< "Cloud name: " << this->name() << nl
987  << " Number of dsmc particles = "
988  << nDSMCParticles
989  << endl;
990 
991  if (nDSMCParticles)
992  {
993  Info<< " Number of molecules = "
994  << nMol << nl
995  << " Mass in system = "
996  << returnReduce(massInSystem(), sumOp<scalar>()) << nl
997  << " Average linear momentum = "
998  << linearMomentum/nMol << nl
999  << " |Average linear momentum| = "
1000  << mag(linearMomentum)/nMol << nl
1001  << " Average linear kinetic energy = "
1002  << linearKineticEnergy/nMol << nl
1003  << " Average internal energy = "
1004  << internalEnergy/nMol << nl
1005  << " Average total energy = "
1006  << (internalEnergy + linearKineticEnergy)/nMol
1007  << endl;
1008  }
1009 }
1010 
1011 
1012 template<class ParcelType>
1015  scalar temperature,
1016  scalar mass
1017 )
1018 {
1019  return
1020  sqrt(physicoChemical::k.value()*temperature/mass)
1021  *rndGen_.sampleNormal<vector>();
1022 }
1023 
1024 
1025 template<class ParcelType>
1028  scalar temperature,
1029  direction iDof
1030 )
1031 {
1032  scalar Ei = 0.0;
1033 
1034  if (iDof == 0)
1035  {
1036  return Ei;
1037  }
1038  else if (iDof == 2)
1039  {
1040  // Special case for iDof = 2, i.e. diatomics;
1041  Ei = -log(rndGen_.scalar01())*physicoChemical::k.value()*temperature;
1042  }
1043  else
1044  {
1045  scalar a = 0.5*iDof - 1;
1046  scalar energyRatio;
1047  scalar P = -1;
1048 
1049  do
1050  {
1051  energyRatio = 10*rndGen_.scalar01();
1052  P = pow((energyRatio/a), a)*exp(a - energyRatio);
1053  } while (P < rndGen_.scalar01());
1054 
1055  Ei = energyRatio*physicoChemical::k.value()*temperature;
1056  }
1057 
1058  return Ei;
1059 }
1060 
1061 
1062 template<class ParcelType>
1064 {
1065  OFstream pObj
1066  (
1067  this->db().time().path()/"parcelPositions_"
1068  + this->name() + "_"
1069  + this->db().time().timeName() + ".obj"
1070  );
1071 
1072  forAllConstIter(typename DSMCCloud<ParcelType>, *this, iter)
1073  {
1074  const ParcelType& p = iter();
1075 
1076  pObj<< "v " << p.position().x()
1077  << " " << p.position().y()
1078  << " " << p.position().z()
1079  << nl;
1080  }
1081 
1082  pObj.flush();
1083 }
1084 
1085 
1086 template<class ParcelType>
1088 {
1090 
1091  // Update the cell occupancy field
1092  cellOccupancy_.setSize(mesh_.nCells());
1093  buildCellOccupancy();
1094 
1095  // Update the inflow BCs
1096  this->inflowBoundary().autoMap(mapper);
1097 }
1098 
1099 
1100 // ************************************************************************* //
Collection of constants.
tetrahedron< point, const point & > tetPointRef
Definition: tetPointRef.H:44
virtual void autoMap(const mapPolyMesh &)
Remap the particles to the correct cells following mesh change.
Definition: DSMCCloud.C:1087
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
tUEqn clear()
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
dimensionedScalar log(const dimensionedScalar &ds)
void readFields(const typename GeoFieldType::Mesh &mesh, const IOobjectList &objects, const HashSet< word > &selectedFields, LIFOStack< regIOobject *> &storedObjects)
Read the selected GeometricFields of the specified type.
Definition: ReadFields.C:244
error FatalError
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:323
void info() const
Print cloud information.
Definition: DSMCCloud.C:970
Templated inflow boundary model class.
Definition: DSMCCloud.H:61
#define forAllIter(Container, container, iter)
Iterate across all elements in the container object of type.
Definition: UList.H:459
Output to file stream.
Definition: OFstream.H:82
uint8_t direction
Definition: direction.H:45
static int myProcNo(const label communicator=0)
Number of this process (starting from masterNo() = 0)
Definition: UPstream.H:429
Templated DSMC particle collision class.
Definition: DSMCCloud.H:55
dimensioned< vector > dimensionedVector
Dimensioned vector obtained from generic dimensioned type.
dimensionedScalar sqrt(const dimensionedScalar &ds)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
const dimensionSet dimless
const dimensionedScalar c
Speed of light in a vacuum.
void autoMap(const mapPolyMesh &)
Remap the cells of particles corresponding to the.
Definition: Cloud.C:341
Class containing mesh-to-mesh mapping information after a change in polyMesh topology.
Definition: mapPolyMesh.H:158
static List< tetIndices > cellTetIndices(const polyMesh &mesh, label cI)
Return the tet decomposition of the given cell, see.
tmp< DimensionedField< TypeR, GeoMesh > > New(const tmp< DimensionedField< TypeR, GeoMesh >> &tdf1, const word &name, const dimensionSet &dimensions)
IOdictionary is derived from dictionary and IOobject to give the dictionary automatic IO functionalit...
Definition: IOdictionary.H:53
Dimension set for the base types.
Definition: dimensionSet.H:120
dynamicFvMesh & mesh
dimensionedScalar exp(const dimensionedScalar &ds)
A class for handling words, derived from string.
Definition: word.H:59
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
bool isType(const Type &t)
Check the typeid.
Definition: typeInfo.H:126
virtual void flush()
Flush stream.
Definition: OSstream.C:207
const word & constant() const
Return constant name.
Definition: TimePaths.H:123
const Type & value() const
Return const reference to value.
word timeName
Definition: getTimeIndex.H:3
static const zero Zero
Definition: zero.H:97
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
Definition: pEqn.H:29
errorManip< error > abort(error &err)
Definition: errorManip.H:131
scalar equipartitionInternalEnergy(scalar temperature, direction internalDegreesOfFreedom)
Generate a random internal energy, sampled from the.
Definition: DSMCCloud.C:1027
label findMax(const ListType &, const label start=0)
Find index of max element (and larger than given element).
dimensionedScalar pos0(const dimensionedScalar &ds)
static const char nl
Definition: Ostream.H:260
const Time & time() const
Return time.
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurrence of given element and return index,.
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
DSMCCloud(const word &cloudName, const fvMesh &mesh, bool readFields=true)
Construct given name and mesh, will read Parcels and fields from.
Definition: DSMCCloud.C:470
void addNewParcel(const vector &position, const label celli, const vector &U, const scalar Ei, const label typeId)
Add new parcel.
Definition: DSMCCloud.C:454
virtual ~DSMCCloud()
Destructor.
Definition: DSMCCloud.C:933
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
void move(TrackCloudType &cloud, typename ParticleType::trackingData &td, const scalar trackTime)
Move the particles.
Definition: Cloud.C:142
vector point
Point is a vector.
Definition: point.H:41
const dimensionedScalar k
Boltzmann constant.
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:78
void evolve()
Evolve the cloud (move, collide)
Definition: DSMCCloud.C:940
messageStream Info
dimensioned< scalar > mag(const dimensioned< Type > &)
Field< vector > vectorField
Specialisation of Field<T> for vector.
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
const word cloudName(propsDict.lookup("cloudName"))
Templated wall interaction model class.
Definition: DSMCCloud.H:58
void dumpParticlePositions() const
Dump particle positions to .obj file.
Definition: DSMCCloud.C:1063
vector equipartitionLinearVelocity(scalar temperature, scalar mass)
Generate a random velocity sampled from the Maxwellian speed.
Definition: DSMCCloud.C:1014
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:92
Templated base class for dsmc cloud.
Definition: DSMCCloud.H:75
fileName path(UMean.rootPath()/UMean.caseName()/functionObjects::writeFile::outputPrefix/"graphs"/UMean.instance())