41 template<
class ParcelType>
44 Info<<
nl <<
"Constructing constant properties for" <<
endl;
45 constProps_.setSize(typeIdList_.size());
47 dictionary moleculeProperties
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"));
95 const dictionary& numberDensitiesDict
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(),
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>();
688 template<
class ParcelType>
693 const IOdictionary& dsmcInitialiseDict
705 mesh_.time().constant(),
707 IOobject::MUST_READ_IF_MODIFIED,
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_",
857 rndGen_(Pstream::myProcNo()),
890 binaryCollisionModel_(),
891 wallInteractionModel_(),
892 inflowBoundaryModel_()
896 initialise(dsmcInitialiseDict);
902 template<
class ParcelType>
909 template<
class ParcelType>
912 typename ParcelType::trackingData td(*
this);
919 this->dumpParticlePositions();
923 this->inflowBoundary().inflow();
929 buildCellOccupancy();
939 template<
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
982 template<
class ParcelType>
991 *rndGen_.GaussNormal<
vector>();
995 template<
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>());
1031 template<
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()
1053 template<
class ParcelType>
1059 cellOccupancy_.setSize(mesh_.nCells());
1060 buildCellOccupancy();
1063 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.
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.
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.
label k
Boltzmann constant.
const dimensionSet dimless
Dimensionless.
tmp< DimensionedField< TypeR, GeoMesh > > New(const tmp< DimensionedField< TypeR, GeoMesh >> &tdf1, const word &name, const dimensionSet &dimensions, const bool initCopy=false)
Global function forwards to reuseTmpDimensionedField::New.
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 INVALID.
dimensionedScalar exp(const dimensionedScalar &ds)
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)
void reduce(const List< UPstream::commsStruct > &comms, T &value, const BinaryOp &bop, const int tag, const label comm)
Reduce inplace (cf. MPI Allreduce) using specified communication schedule.
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.
void readFields(const typename GeoFieldType::Mesh &mesh, const IOobjectList &objects, const wordHashSet &selectedFields, LIFOStack< regIOobject *> &storedObjects)
Read the selected GeometricFields of the templated type.
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.
static constexpr const zero Zero
Global zero (0)