41template<
class ParcelType>
44 Info<<
nl <<
"Constructing constant properties for" <<
endl;
45 constProps_.setSize(typeIdList_.size());
49 particleProperties_.subDict(
"moleculeProperties")
54 const word&
id = typeIdList_[i];
58 const dictionary& molDict(moleculeProperties.subDict(
id));
60 constProps_[i] =
typename ParcelType::constantProperties(molDict);
65template<
class ParcelType>
68 for (
auto& list : cellOccupancy_)
73 for (ParcelType&
p : *
this)
75 cellOccupancy_[
p.cell()].append(&
p);
80template<
class ParcelType>
88 const scalar temperature
90 dsmcInitialiseDict.
get<scalar>(
"temperature")
97 dsmcInitialiseDict.
subDict(
"numberDensities")
100 List<word> molecules(numberDensitiesDict.toc());
106 numberDensities[i] = numberDensitiesDict.get<scalar>(molecules[i]);
109 numberDensities /= nParticle_;
111 forAll(mesh_.cells(), celli)
123 scalar tetVolume = tet.
mag();
127 const word& moleculeName(molecules[i]);
129 label typeId = typeIdList_.find(moleculeName);
134 <<
"typeId " << moleculeName <<
"not defined." <<
nl
135 <<
abort(FatalError);
138 const typename ParcelType::constantProperties& cP =
141 scalar numberDensity = numberDensities[i];
144 scalar particlesRequired = numberDensity*tetVolume;
147 label nParticlesToInsert = label(particlesRequired);
153 (particlesRequired - nParticlesToInsert)
154 > rndGen_.sample01<scalar>()
157 nParticlesToInsert++;
160 for (label pI = 0; pI < nParticlesToInsert; pI++)
164 vector U = equipartitionLinearVelocity
170 scalar Ei = equipartitionInternalEnergy
173 cP.internalDegreesOfFreedom()
178 addNewParcel(
p, celli,
U, Ei, typeId);
188 label mostAbundantType(
findMax(numberDensities));
190 const typename ParcelType::constantProperties& cP = constProps
195 sigmaTcRMax_.primitiveFieldRef() = cP.sigmaT()*maxwellianMostProbableSpeed
201 sigmaTcRMax_.correctBoundaryConditions();
205template<
class ParcelType>
208 if (!binaryCollision().active())
216 scalar deltaT =
mesh().time().deltaTValue();
218 label collisionCandidates = 0;
220 label collisions = 0;
222 forAll(cellOccupancy_, celli)
226 label nC(cellParcels.size());
242 const point& cC = mesh_.cellCentres()[celli];
246 const ParcelType&
p = *cellParcels[i];
247 vector relPos =
p.position() - cC;
250 pos0(relPos.x()) + 2*
pos0(relPos.y()) + 4*
pos0(relPos.z());
252 subCells[subCell].append(i);
253 whichSubCell[i] = subCell;
258 scalar sigmaTcRMax = sigmaTcRMax_[celli];
260 scalar selectedPairs =
261 collisionSelectionRemainder_[celli]
262 + 0.5*nC*(nC - 1)*nParticle_*sigmaTcRMax*deltaT
263 /mesh_.cellVolumes()[celli];
265 label nCandidates(selectedPairs);
266 collisionSelectionRemainder_[celli] = selectedPairs - nCandidates;
267 collisionCandidates += nCandidates;
269 for (label c = 0;
c < nCandidates;
c++)
275 label candidateP = rndGen_.position<label>(0, nC - 1);
278 label candidateQ = -1;
280 List<label> subCellPs = subCells[whichSubCell[candidateP]];
281 label nSC = subCellPs.
size();
291 label i = rndGen_.position<label>(0, nSC - 1);
292 candidateQ = subCellPs[i];
293 }
while (candidateP == candidateQ);
303 candidateQ = rndGen_.position<label>(0, nC - 1);
304 }
while (candidateP == candidateQ);
324 ParcelType& parcelP = *cellParcels[candidateP];
325 ParcelType& parcelQ = *cellParcels[candidateQ];
327 scalar sigmaTcR = binaryCollision().sigmaTcR
337 if (sigmaTcR > sigmaTcRMax_[celli])
339 sigmaTcRMax_[celli] = sigmaTcR;
342 if ((sigmaTcR/sigmaTcRMax) > rndGen_.sample01<scalar>())
344 binaryCollision().collide
360 sigmaTcRMax_.correctBoundaryConditions();
362 if (collisionCandidates)
364 Info<<
" Collisions = "
366 <<
" Acceptance rate = "
367 << scalar(collisions)/scalar(collisionCandidates) <<
nl
377template<
class ParcelType>
393template<
class ParcelType>
398 scalarField& dsmcRhoN = dsmcRhoN_.primitiveFieldRef();
399 scalarField& linearKE = linearKE_.primitiveFieldRef();
400 scalarField& internalE = internalE_.primitiveFieldRef();
402 vectorField& momentum = momentum_.primitiveFieldRef();
404 for (
const ParcelType&
p : *
this)
406 const label celli =
p.cell();
409 rhoM[celli] += constProps(
p.typeId()).mass();
411 linearKE[celli] += 0.5*constProps(
p.typeId()).mass()*(
p.U() &
p.U());
412 internalE[celli] +=
p.Ei();
413 iDof[celli] += constProps(
p.typeId()).internalDegreesOfFreedom();
414 momentum[celli] += constProps(
p.typeId()).mass()*
p.U();
417 rhoN *= nParticle_/
mesh().cellVolumes();
418 rhoN_.correctBoundaryConditions();
420 rhoM *= nParticle_/
mesh().cellVolumes();
421 rhoM_.correctBoundaryConditions();
423 dsmcRhoN_.correctBoundaryConditions();
425 linearKE *= nParticle_/
mesh().cellVolumes();
426 linearKE_.correctBoundaryConditions();
428 internalE *= nParticle_/
mesh().cellVolumes();
429 internalE_.correctBoundaryConditions();
431 iDof *= nParticle_/
mesh().cellVolumes();
432 iDof_.correctBoundaryConditions();
434 momentum *= nParticle_/
mesh().cellVolumes();
435 momentum_.correctBoundaryConditions();
441template<
class ParcelType>
451 this->addParticle(
new ParcelType(mesh_, position, celli,
U, Ei, typeId));
457template<
class ParcelType>
480 typeIdList_(particleProperties_.
lookup(
"typeIdList")),
481 nParticle_(particleProperties_.get<scalar>(
"nEquivalentParticles")),
482 cellOccupancy_(mesh_.nCells()),
487 this->
name() +
"SigmaTcRMax",
495 collisionSelectionRemainder_
499 this->
name() +
":collisionSelectionRemainder",
646 binaryCollisionModel_
654 wallInteractionModel_
672 buildCellOccupancy();
676 forAll(collisionSelectionRemainder_, i)
678 collisionSelectionRemainder_[i] = rndGen_.
sample01<scalar>();
688template<
class ParcelType>
711 typeIdList_(particleProperties_.
lookup(
"typeIdList")),
712 nParticle_(particleProperties_.get<scalar>(
"nEquivalentParticles")),
718 this->
name() +
"SigmaTcRMax",
726 zeroGradientFvPatchScalarField::typeName
728 collisionSelectionRemainder_
732 this->
name() +
":collisionSelectionRemainder",
756 this->
name() +
"fD_",
769 this->
name() +
"rhoN_",
782 this->
name() +
"rhoM_",
795 this->
name() +
"dsmcRhoN_",
808 this->
name() +
"linearKE_",
821 this->
name() +
"internalE_",
834 this->
name() +
"iDof_",
847 this->
name() +
"momentum_",
890 binaryCollisionModel_(),
891 wallInteractionModel_(),
892 inflowBoundaryModel_()
896 initialise(dsmcInitialiseDict);
902template<
class ParcelType>
909template<
class ParcelType>
919 this->dumpParticlePositions();
923 this->inflowBoundary().inflow();
929 buildCellOccupancy();
939template<
class ParcelType>
942 label nDSMCParticles = this->size();
945 scalar nMol = nDSMCParticles*nParticle_;
947 vector linearMomentum = linearMomentumOfSystem();
950 scalar linearKineticEnergy = linearKineticEnergyOfSystem();
953 scalar internalEnergy = internalEnergyOfSystem();
957 <<
" Number of dsmc particles = "
963 Info<<
" Number of molecules = "
965 <<
" Mass in system = "
967 <<
" Average linear momentum = "
968 << linearMomentum/nMol <<
nl
969 <<
" |Average linear momentum| = "
970 <<
mag(linearMomentum)/nMol <<
nl
971 <<
" Average linear kinetic energy = "
972 << linearKineticEnergy/nMol <<
nl
973 <<
" Average internal energy = "
974 << internalEnergy/nMol <<
nl
975 <<
" Average total energy = "
976 << (internalEnergy + linearKineticEnergy)/nMol
982template<
class ParcelType>
991 *rndGen_.GaussNormal<
vector>();
995template<
class ParcelType>
1011 -
log(rndGen_.sample01<scalar>())
1017 const scalar a = 0.5*iDof - 1;
1018 scalar energyRatio = 0;
1023 energyRatio = 10*rndGen_.sample01<scalar>();
1024 P =
pow((energyRatio/a), a)*
exp(a - energyRatio);
1025 }
while (P < rndGen_.sample01<scalar>());
1031template<
class ParcelType>
1036 this->db().time().
path()/
"parcelPositions_"
1037 + this->
name() +
"_"
1038 + this->db().time().
timeName() +
".obj"
1041 for (
const ParcelType&
p : *
this)
1043 pObj<<
"v " <<
p.position().x()
1044 <<
' ' <<
p.position().y()
1045 <<
' ' <<
p.position().z()
1053template<
class ParcelType>
1059 cellOccupancy_.setSize(mesh_.nCells());
1060 buildCellOccupancy();
1063 this->inflowBoundary().autoMap(mapper);
reduce(hasMovingMesh, orOp< bool >())
Templated DSMC particle collision class.
Base cloud calls templated on particle type.
void autoMap(const mapPolyMesh &)
Remap the cells of particles corresponding to the.
Virtual abstract base class for templated DSMCCloud.
Templated base class for dsmc cloud.
virtual ~DSMCCloud()
Destructor.
vector equipartitionLinearVelocity(scalar temperature, scalar mass)
Generate a random velocity sampled from the Maxwellian speed.
virtual void autoMap(const mapPolyMesh &)
Remap the particles to the correct cells following mesh change.
void evolve()
Evolve the cloud (move, collide)
scalar equipartitionInternalEnergy(scalar temperature, direction internalDegreesOfFreedom)
Generate a random internal energy, sampled from the.
void clear()
Clear the Cloud.
void addNewParcel(const vector &position, const label celli, const vector &U, const scalar Ei, const label typeId)
Add new parcel.
void dumpParticlePositions() const
Dump particle positions to .obj file.
void info() const
Print cloud information.
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects.
Generic templated field type.
IOdictionary is derived from dictionary and IOobject to give the dictionary automatic IO functionalit...
Defines the attributes of an object for which implicit objectRegistry management is supported,...
Templated inflow boundary model class.
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
Output to file stream, using an OSstream.
virtual void flush()
Flush stream.
Inter-processor communications stream.
Type sample01()
Return a sample whose components lie in the range [0,1].
void size(const label n)
Older name for setAddressableSize.
Templated wall interaction model class.
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
T get(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
const dictionary & subDict(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
Find and return a sub-dictionary.
Dimension set for the base types, which can be used to implement rigorous dimension checking for alge...
const Type & value() const
Return const reference to value.
Class used to pass data into container.
Reads fields from the time directories and adds them to the mesh database for further post-processing...
Mesh data needed to do the Finite Volume discretisation.
Class containing mesh-to-mesh mapping information after a change in polyMesh topology.
constant condensation/saturation model.
Lookup type of boundary radiation properties.
Storage and named access for the indices of a tet which is part of the decomposition of a cell.
tetPointRef tet(const polyMesh &mesh) const
Return the geometry corresponding to this tet.
scalar mag() const
Return volume.
Point randomPoint(Random &rndGen) const
Return a random point in the tetrahedron from a.
void initialise()
Initialise integral-scale box properties.
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
A class for handling words, derived from Foam::string.
fileName path(UMean.rootPath()/UMean.caseName()/"graphs"/UMean.instance())
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
const dimensionedScalar k
Boltzmann constant.
const dimensionedScalar c
Speed of light in a vacuum.
Different types of constants.
dimensionedScalar exp(const dimensionedScalar &ds)
const dimensionSet dimless
Dimensionless.
dimensionedScalar pos0(const dimensionedScalar &ds)
messageStream Info
Information stream (stdout output on master, null elsewhere)
void readFields(const typename GeoFieldType::Mesh &mesh, const IOobjectList &objects, const wordHashSet &selectedFields, LIFOStack< regIOobject * > &storedObjects)
Read the selected GeometricFields of the templated type.
dimensionedScalar log(const dimensionedScalar &ds)
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
Ostream & endl(Ostream &os)
Add newline and flush stream.
dimensionedScalar sqrt(const dimensionedScalar &ds)
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
void reduce(const List< UPstream::commsStruct > &comms, T &value, const BinaryOp &bop, const int tag, const label comm)
errorManip< error > abort(error &err)
static constexpr const zero Zero
Global zero (0)
tmp< DimensionedField< TypeR, GeoMesh > > New(const tmp< DimensionedField< TypeR, GeoMesh > > &tdf1, const word &name, const dimensionSet &dimensions)
Global function forwards to reuseTmpDimensionedField::New.
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
T returnReduce(const T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Reduce (copy) and return value.
label findMax(const ListType &input, label start=0)
constexpr char nl
The newline '\n' character (0x0a)
#define forAll(list, i)
Loop across all elements in list.
const word cloudName(propsDict.get< word >("cloud"))