41 template<
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);
65 template<
class ParcelType>
68 for (
auto& list : cellOccupancy_)
73 for (ParcelType&
p : *
this)
75 cellOccupancy_[
p.cell()].append(&
p);
80 template<
class ParcelType>
83 const IOdictionary& dsmcInitialiseDict
88 const scalar temperature
90 dsmcInitialiseDict.get<scalar>(
"temperature")
93 const vector velocity(dsmcInitialiseDict.get<
vector>(
"velocity"));
97 dsmcInitialiseDict.subDict(
"numberDensities")
100 List<word> molecules(numberDensitiesDict.toc());
102 Field<scalar> numberDensities(molecules.size());
106 numberDensities[i] = numberDensitiesDict.get<scalar>(molecules[i]);
109 numberDensities /= nParticle_;
111 forAll(mesh_.cells(), celli)
121 const tetIndices& cellTetIs = cellTets[tetI];
123 scalar tetVolume = tet.mag();
127 const word& moleculeName(molecules[i]);
129 label typeId = typeIdList_.find(moleculeName);
134 <<
"typeId " << moleculeName <<
"not defined." <<
nl 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++)
162 point p = tet.randomPoint(rndGen_);
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();
205 template<
class ParcelType>
208 if (!binaryCollision().active())
214 List<DynamicList<label>> subCells(8);
216 scalar deltaT =
mesh().time().deltaTValue();
218 label collisionCandidates = 0;
220 label collisions = 0;
222 forAll(cellOccupancy_, celli)
224 const DynamicList<ParcelType*>& cellParcels(cellOccupancy_[celli]);
226 label nC(cellParcels.size());
240 List<label> whichSubCell(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
356 reduce(collisions, sumOp<label>());
358 reduce(collisionCandidates, sumOp<label>());
360 sigmaTcRMax_.correctBoundaryConditions();
362 if (collisionCandidates)
364 Info<<
" Collisions = " 366 <<
" Acceptance rate = " 367 << scalar(collisions)/scalar(collisionCandidates) <<
nl 377 template<
class ParcelType>
393 template<
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();
441 template<
class ParcelType>
451 this->addParticle(
new ParcelType(mesh_, position, celli,
U, Ei, typeId));
457 template<
class ParcelType>
474 mesh_.time().constant(),
481 typeIdList_(particleProperties_.
lookup(
"typeIdList")),
482 nParticle_(particleProperties_.
get<scalar>(
"nEquivalentParticles")),
483 cellOccupancy_(mesh_.nCells()),
488 this->
name() +
"SigmaTcRMax",
497 collisionSelectionRemainder_
501 IOobject::scopedName(this->
name(),
"collisionSelectionRemainder"),
653 binaryCollisionModel_
661 wallInteractionModel_
679 buildCellOccupancy();
683 forAll(collisionSelectionRemainder_, i)
685 collisionSelectionRemainder_[i] = rndGen_.
sample01<scalar>();
695 template<
class ParcelType>
700 const IOdictionary& dsmcInitialiseDict
712 mesh_.time().constant(),
714 IOobject::READ_MODIFIED,
719 typeIdList_(particleProperties_.lookup(
"typeIdList")),
720 nParticle_(particleProperties_.
get<scalar>(
"nEquivalentParticles")),
726 this->
name() +
"SigmaTcRMax",
730 IOobject::AUTO_WRITE,
737 collisionSelectionRemainder_
741 IOobject::scopedName(this->
name(),
"collisionSelectionRemainder"),
765 this->
name() +
"fD_",
778 this->
name() +
"rhoN_",
791 this->
name() +
"rhoM_",
804 this->
name() +
"dsmcRhoN_",
817 this->
name() +
"linearKE_",
830 this->
name() +
"internalE_",
843 this->
name() +
"iDof_",
856 this->
name() +
"momentum_",
866 rndGen_(Pstream::myProcNo()),
899 binaryCollisionModel_(),
900 wallInteractionModel_(),
901 inflowBoundaryModel_()
905 initialise(dsmcInitialiseDict);
911 template<
class ParcelType>
918 template<
class ParcelType>
921 typename ParcelType::trackingData
td(*
this);
928 this->dumpParticlePositions();
932 this->inflowBoundary().inflow();
938 buildCellOccupancy();
948 template<
class ParcelType>
951 label nDSMCParticles = this->size();
954 scalar nMol = nDSMCParticles*nParticle_;
956 vector linearMomentum = linearMomentumOfSystem();
959 scalar linearKineticEnergy = linearKineticEnergyOfSystem();
962 scalar internalEnergy = internalEnergyOfSystem();
966 <<
" Number of dsmc particles = " 972 Info<<
" Number of molecules = " 974 <<
" Mass in system = " 976 <<
" Average linear momentum = " 977 << linearMomentum/nMol <<
nl 978 <<
" |Average linear momentum| = " 979 <<
mag(linearMomentum)/nMol <<
nl 980 <<
" Average linear kinetic energy = " 981 << linearKineticEnergy/nMol <<
nl 982 <<
" Average internal energy = " 983 << internalEnergy/nMol <<
nl 984 <<
" Average total energy = " 985 << (internalEnergy + linearKineticEnergy)/nMol
991 template<
class ParcelType>
1000 *rndGen_.GaussNormal<
vector>();
1004 template<
class ParcelType>
1020 -
log(rndGen_.sample01<scalar>())
1026 const scalar a = 0.5*iDof - 1;
1027 scalar energyRatio = 0;
1032 energyRatio = 10*rndGen_.sample01<scalar>();
1033 P =
pow((energyRatio/a), a)*
exp(a - energyRatio);
1034 }
while (P < rndGen_.sample01<scalar>());
1040 template<
class ParcelType>
1045 this->db().time().
path()/
"parcelPositions_" 1046 + this->
name() +
"_" 1047 + this->db().time().
timeName() +
".obj" 1050 for (
const ParcelType&
p : *
this)
1052 pObj<<
"v " <<
p.position().x()
1053 <<
' ' <<
p.position().y()
1054 <<
' ' <<
p.position().z()
1062 template<
class ParcelType>
1068 cellOccupancy_.setSize(mesh_.nCells());
1069 buildCellOccupancy();
1072 this->inflowBoundary().autoMap(mapper);
Different types of constants.
List< ReturnType > get(const UPtrList< T > &list, const AccessOp &aop)
List of values generated by applying the access operation to each list item.
label findMax(const ListType &input, label start=0)
Linear search for the index of the max element, similar to std::max_element but for lists and returns...
virtual void autoMap(const mapPolyMesh &)
Remap the particles to the correct cells following mesh change.
const word zeroGradientType
A zeroGradient patch field type.
dimensionedScalar log(const dimensionedScalar &ds)
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
error FatalError
Error stream (stdout output on all processes), with additional 'FOAM FATAL ERROR' header text and sta...
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
void info() const
Print cloud information.
tetrahedron< point, const point & > tetPointRef
A tetrahedron using referred points.
wallPoints::trackData td(isBlockedFace, regionToBlockSize)
Templated inflow boundary model class.
constexpr char nl
The newline '\n' character (0x0a)
Templated DSMC particle collision class.
dimensioned< vector > dimensionedVector
Dimensioned vector obtained from generic dimensioned type.
dimensionedScalar sqrt(const dimensionedScalar &ds)
Ostream & endl(Ostream &os)
Add newline and flush stream.
tmp< DimensionedField< TypeR, GeoMesh > > New(const tmp< DimensionedField< TypeR, GeoMesh >> &tf1, const word &name, const dimensionSet &dimensions, const bool initCopy=false)
Global function forwards to reuseTmpDimensionedField::New.
label k
Boltzmann constant.
const dimensionSet dimless
Dimensionless.
Lookup type of boundary radiation properties.
GeometricField< vector, fvPatchField, volMesh > volVectorField
T returnReduce(const T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Perform reduction on a copy, using specified binary operation.
void autoMap(const mapPolyMesh &)
Remap the cells of particles corresponding to the.
static List< tetIndices > cellTetIndices(const polyMesh &mesh, label cI)
Return the tet decomposition of the given cell, see.
#define forAll(list, i)
Loop across all elements in list.
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Type sample01()
Return a sample whose components lie in the range [0,1].
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for expressions::valueTypeCode::INVALID.
dimensionedScalar exp(const dimensionedScalar &ds)
void readFields(const typename GeoFieldType::Mesh &mesh, const IOobjectList &objects, const NameMatchPredicate &selectedFields, DynamicList< regIOobject *> &storedObjects)
Read the selected GeometricFields of the templated type and store on the objectRegistry.
const word cloudName(propsDict.get< word >("cloud"))
A class for handling words, derived from Foam::string.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Inter-processor communications stream.
errorManip< error > abort(error &err)
Base cloud calls templated on particle type.
scalar equipartitionInternalEnergy(scalar temperature, direction internalDegreesOfFreedom)
Generate a random internal energy, sampled from the.
void clear()
Clear the Cloud.
dimensionedScalar pos0(const dimensionedScalar &ds)
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
int debug
Static debugging option.
fileName path(UMean.rootPath()/UMean.caseName()/"graphs"/UMean.instance())
void addNewParcel(const vector &position, const label celli, const vector &U, const scalar Ei, const label typeId)
Add new parcel.
virtual ~DSMCCloud()
Destructor.
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
void move(TrackCloudType &cloud, typename ParticleType::trackingData &td, const scalar trackTime)
Move the particles.
vector point
Point is a vector.
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
Mesh data needed to do the Finite Volume discretisation.
const dimensionedScalar c
Speed of light in a vacuum.
void evolve()
Evolve the cloud (move, collide)
messageStream Info
Information stream (stdout output on master, null elsewhere)
Field< vector > vectorField
Specialisation of Field<T> for vector.
Virtual abstract base class for templated DSMCCloud.
Templated wall interaction model class.
void dumpParticlePositions() const
Dump particle positions to .obj file.
vector equipartitionLinearVelocity(scalar temperature, scalar mass)
Generate a random velocity sampled from the Maxwellian speed.
Defines the attributes of an object for which implicit objectRegistry management is supported...
Templated base class for dsmc cloud.
void reduce(T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Reduce inplace (cf. MPI Allreduce) using linear/tree communication schedule.
static constexpr const zero Zero
Global zero (0)