38 template<
class ParcelType>
41 Info<<
nl <<
"Constructing constant properties for" <<
endl;
42 constProps_.setSize(typeIdList_.size());
44 dictionary moleculeProperties
46 particleProperties_.subDict(
"moleculeProperties")
51 const word& id(typeIdList_[i]);
55 const dictionary& molDict(moleculeProperties.subDict(
id));
58 typename ParcelType::constantProperties(molDict);
63 template<
class ParcelType>
68 cellOccupancy_[cO].clear();
71 forAllIter(
typename DSMCCloud<ParcelType>, *
this, iter)
73 cellOccupancy_[iter().cell()].append(&iter());
78 template<
class ParcelType>
81 const IOdictionary& dsmcInitialiseDict
86 const scalar temperature
88 dsmcInitialiseDict.template lookup<scalar>(
"temperature")
91 const vector velocity(dsmcInitialiseDict.lookup(
"velocity"));
93 const dictionary& numberDensitiesDict
95 dsmcInitialiseDict.subDict(
"numberDensities")
98 List<word> molecules(numberDensitiesDict.toc());
100 Field<scalar> numberDensities(molecules.size());
105 numberDensitiesDict.lookup<scalar>(molecules[i]);
108 numberDensities /= nParticle_;
110 forAll(mesh_.cells(), celli)
120 const tetIndices& cellTetIs = cellTets[tetI];
122 scalar tetVolume = tet.mag();
126 const word& moleculeName(molecules[i]);
133 <<
"typeId " << moleculeName <<
"not defined." <<
nl
137 const typename ParcelType::constantProperties& cP =
140 scalar numberDensity = numberDensities[i];
143 scalar particlesRequired = numberDensity*tetVolume;
146 label nParticlesToInsert =
label(particlesRequired);
152 (particlesRequired - nParticlesToInsert)
156 nParticlesToInsert++;
159 for (
label pI = 0; pI < nParticlesToInsert; pI++)
161 point p = tet.randomPoint(rndGen_);
163 vector U = equipartitionLinearVelocity
169 scalar Ei = equipartitionInternalEnergy
172 cP.internalDegreesOfFreedom()
177 addNewParcel(
p, celli,
U, Ei, typeId);
189 const typename ParcelType::constantProperties& cP = constProps
194 sigmaTcRMax_.primitiveFieldRef() = cP.sigmaT()*maxwellianMostProbableSpeed
200 sigmaTcRMax_.correctBoundaryConditions();
204 template<
class ParcelType>
207 if (
isType<NoBinaryCollision<DSMCCloud<ParcelType>>>(binaryCollision()))
213 List<DynamicList<label>> subCells(8);
215 scalar deltaT = mesh().time().deltaTValue();
217 label collisionCandidates = 0;
219 label collisions = 0;
221 forAll(cellOccupancy_, celli)
223 const DynamicList<ParcelType*>& cellParcels(cellOccupancy_[celli]);
225 label nC(cellParcels.size());
239 List<label> whichSubCell(cellParcels.size());
241 const point& cC = mesh_.cellCentres()[celli];
245 const ParcelType&
p = *cellParcels[i];
246 vector relPos =
p.position(mesh()) - cC;
249 pos0(relPos.x()) + 2*
pos0(relPos.y()) + 4*
pos0(relPos.z());
251 subCells[subCell].append(i);
252 whichSubCell[i] = subCell;
257 scalar sigmaTcRMax = sigmaTcRMax_[celli];
259 scalar selectedPairs =
260 collisionSelectionRemainder_[celli]
261 + 0.5*nC*(nC - 1)*nParticle_*sigmaTcRMax*deltaT
262 /mesh_.cellVolumes()[celli];
264 label nCandidates(selectedPairs);
265 collisionSelectionRemainder_[celli] = selectedPairs - nCandidates;
266 collisionCandidates += nCandidates;
268 for (
label c = 0;
c < nCandidates;
c++)
274 label candidateP = rndGen_.sampleAB<
label>(0, nC);
277 label candidateQ = -1;
279 List<label> subCellPs = subCells[whichSubCell[candidateP]];
280 label nSC = subCellPs.size();
290 candidateQ = subCellPs[rndGen_.sampleAB<
label>(0, nSC)];
291 }
while (candidateP == candidateQ);
301 candidateQ = rndGen_.sampleAB<
label>(0, nC);
302 }
while (candidateP == candidateQ);
322 ParcelType& parcelP = *cellParcels[candidateP];
323 ParcelType& parcelQ = *cellParcels[candidateQ];
325 scalar sigmaTcR = binaryCollision().sigmaTcR
335 if (sigmaTcR > sigmaTcRMax_[celli])
337 sigmaTcRMax_[celli] = sigmaTcR;
340 if ((sigmaTcR/sigmaTcRMax) > rndGen_.scalar01())
342 binaryCollision().collide
354 reduce(collisions, sumOp<label>());
356 reduce(collisionCandidates, sumOp<label>());
358 sigmaTcRMax_.correctBoundaryConditions();
360 if (collisionCandidates)
362 Info<<
" Collisions = "
364 <<
" Acceptance rate = "
365 << scalar(collisions)/scalar(collisionCandidates) <<
nl
375 template<
class ParcelType>
383 dimensionSet(1, -1, -2, 0, 0),
397 dimensionSet(1, -2, -1, 0, 0),
403 template<
class ParcelType>
408 scalarField& dsmcRhoN = dsmcRhoN_.primitiveFieldRef();
409 scalarField& linearKE = linearKE_.primitiveFieldRef();
410 scalarField& internalE = internalE_.primitiveFieldRef();
412 vectorField& momentum = momentum_.primitiveFieldRef();
416 const ParcelType&
p = iter();
417 const label celli =
p.cell();
420 rhoM[celli] += constProps(
p.typeId()).mass();
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();
428 rhoN *= nParticle_/mesh().cellVolumes();
429 rhoN_.correctBoundaryConditions();
431 rhoM *= nParticle_/mesh().cellVolumes();
432 rhoM_.correctBoundaryConditions();
434 dsmcRhoN_.correctBoundaryConditions();
436 linearKE *= nParticle_/mesh().cellVolumes();
437 linearKE_.correctBoundaryConditions();
439 internalE *= nParticle_/mesh().cellVolumes();
440 internalE_.correctBoundaryConditions();
442 iDof *= nParticle_/mesh().cellVolumes();
443 iDof_.correctBoundaryConditions();
445 momentum *= nParticle_/mesh().cellVolumes();
446 momentum_.correctBoundaryConditions();
452 template<
class ParcelType>
462 this->addParticle(
new ParcelType(mesh_, position, celli,
U, Ei, typeId));
468 template<
class ParcelType>
484 mesh_.time().constant(),
490 typeIdList_(particleProperties_.lookup(
"typeIdList")),
493 particleProperties_.template lookup<scalar>(
"nEquivalentParticles")
495 cellOccupancy_(mesh_.nCells()),
500 this->
name() +
"SigmaTcRMax",
508 collisionSelectionRemainder_
512 this->
name() +
":collisionSelectionRemainder",
659 binaryCollisionModel_
667 wallInteractionModel_
685 buildCellOccupancy();
689 forAll(collisionSelectionRemainder_, i)
691 collisionSelectionRemainder_[i] = rndGen_.
scalar01();
701 template<
class ParcelType>
717 mesh_.time().constant(),
723 typeIdList_(particleProperties_.lookup(
"typeIdList")),
726 particleProperties_.template lookup<scalar>(
"nEquivalentParticles")
733 this->
name() +
"SigmaTcRMax",
741 zeroGradientFvPatchScalarField::typeName
743 collisionSelectionRemainder_
747 this->
name() +
":collisionSelectionRemainder",
771 this->
name() +
"fD_",
789 this->
name() +
"rhoN_",
802 this->
name() +
"rhoM_",
815 this->
name() +
"dsmcRhoN_",
828 this->
name() +
"linearKE_",
841 this->
name() +
"internalE_",
854 this->
name() +
"iDof_",
867 this->
name() +
"momentum_",
920 binaryCollisionModel_(),
921 wallInteractionModel_(),
922 inflowBoundaryModel_()
926 initialise(dsmcInitialiseDict);
932 template<
class ParcelType>
939 template<
class ParcelType>
942 typename ParcelType::trackingData td(*
this);
949 this->dumpParticlePositions();
953 this->inflowBoundary().inflow();
959 buildCellOccupancy();
969 template<
class ParcelType>
972 label nDSMCParticles = this->size();
975 scalar nMol = nDSMCParticles*nParticle_;
977 vector linearMomentum = linearMomentumOfSystem();
980 scalar linearKineticEnergy = linearKineticEnergyOfSystem();
983 scalar internalEnergy = internalEnergyOfSystem();
987 <<
" Number of dsmc particles = "
993 Info<<
" Number of molecules = "
995 <<
" Mass in system = "
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
1012 template<
class ParcelType>
1021 *rndGen_.sampleNormal<
vector>();
1025 template<
class ParcelType>
1045 scalar a = 0.5*iDof - 1;
1051 energyRatio = 10*rndGen_.scalar01();
1052 P =
pow((energyRatio/a), a)*
exp(a - energyRatio);
1053 }
while (P < rndGen_.scalar01());
1062 template<
class ParcelType>
1067 this->db().time().path()/
"parcelPositions_"
1068 + this->
name() +
"_"
1069 + this->db().time().
name() +
".obj"
1074 const point pos = iter().position(mesh());
1076 pObj<<
"v " <<
pos.x() <<
" " <<
pos.y() <<
" " <<
pos.z() <<
nl;
1083 template<
class ParcelType>
1089 cellOccupancy_.setSize(mesh_.nCells());
1090 buildCellOccupancy();
1093 this->inflowBoundary().topoChange();
1097 template<
class ParcelType>
1103 cellOccupancy_.setSize(mesh_.nCells());
1104 buildCellOccupancy();
1107 this->inflowBoundary().topoChange();
1111 template<
class ParcelType>
1117 cellOccupancy_.setSize(mesh_.nCells());
1118 buildCellOccupancy();
1121 this->inflowBoundary().topoChange();
#define forAll(list, i)
Loop across all elements in list.
#define forAllIter(Container, container, iter)
Iterate across all elements in the container object of type.
#define forAllConstIter(Container, container, iter)
Iterate across all elements in the container object of type.
Templated DSMC particle collision class.
Base cloud calls templated on particle type.
virtual void topoChange(const polyTopoChangeMap &)
Update topology using the given map.
virtual void distribute(const polyDistributionMap &)
Redistribute or update using the given distribution map.
virtual void mapMesh(const polyMeshMap &)
Update from another mesh using the given map.
void move(TrackCloudType &cloud, typename ParticleType::trackingData &td)
Move the particles.
Templated base class for dsmc cloud.
DSMCCloud(const word &cloudName, const fvMesh &mesh, bool readFields=true)
Construct given name and mesh, will read Parcels and fields from.
virtual ~DSMCCloud()
Destructor.
vector equipartitionLinearVelocity(scalar temperature, scalar mass)
Generate a random velocity sampled from the Maxwellian speed.
virtual void topoChange(const polyTopoChangeMap &)
Update topology using the given map.
virtual void distribute(const polyDistributionMap &)
Redistribute or update using the given distribution map.
virtual void mapMesh(const polyMeshMap &)
Update from another mesh using the given map.
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.
Generic GeometricField class.
IOdictionary is derived from dictionary and IOobject to give the dictionary automatic IO functionalit...
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Templated inflow boundary model class.
virtual void flush()
Flush stream.
Inter-processor communications stream.
scalar scalar01()
Advance the state and return a scalar sample from a uniform.
Templated wall interaction model class.
Dimension set for the base types.
Mesh data needed to do the Finite Volume discretisation.
Class containing mesh-to-mesh mapping information after a mesh distribution where we send parts of me...
Class containing mesh-to-mesh mapping information.
static List< tetIndices > cellTetIndices(const polyMesh &mesh, label cI)
Return the tet decomposition of the given cell, see.
Class containing mesh-to-mesh mapping information after a change in polyMesh topology.
A class for handling words, derived from string.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
const dimensionedScalar c
Speed of light in a vacuum.
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
dimensionedScalar pos(const dimensionedScalar &ds)
dimensionedScalar exp(const dimensionedScalar &ds)
dimensionedScalar pos0(const dimensionedScalar &ds)
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
label findMax(const ListType &, const label start=0)
Find index of max element (and larger than given element).
Ostream & endl(Ostream &os)
Add newline and flush stream.
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.
errorManip< error > abort(error &err)
const dimensionSet dimless
bool isType(const Type &t)
Check the typeid.
vector point
Point is a vector.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
dimensionedScalar log(const dimensionedScalar &ds)
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
Vector< scalar > vector
A scalar version of the templated Vector.
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
dimensionedScalar sqrt(const dimensionedScalar &ds)
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
dimensioned< scalar > mag(const dimensioned< Type > &)
tmp< DimensionedField< TypeR, GeoMesh > > New(const tmp< DimensionedField< TypeR, GeoMesh >> &tdf1, const word &name, const dimensionSet &dimensions)
dimensioned< vector > dimensionedVector
Dimensioned vector obtained from generic dimensioned type.
Field< vector > vectorField
Specialisation of Field<T> for vector.
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurrence of given element and return index,.
word name(const complex &)
Return a string representation of a complex.
tetrahedron< point, const point & > tetPointRef
const word cloudName(propsDict.lookup("cloudName"))