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",
640 binaryCollisionModel_
648 wallInteractionModel_
666 buildCellOccupancy();
670 forAll(collisionSelectionRemainder_, i)
672 collisionSelectionRemainder_[i] = rndGen_.
sample01<scalar>();
682 template<
class ParcelType>
687 const IOdictionary& dsmcInitialiseDict
699 mesh_.time().constant(),
701 IOobject::MUST_READ_IF_MODIFIED,
705 typeIdList_(particleProperties_.lookup(
"typeIdList")),
706 nParticle_(particleProperties_.
get<scalar>(
"nEquivalentParticles")),
712 this->
name() +
"SigmaTcRMax",
722 collisionSelectionRemainder_
726 this->
name() +
":collisionSelectionRemainder",
750 this->
name() +
"fD_",
763 this->
name() +
"rhoN_",
776 this->
name() +
"rhoM_",
789 this->
name() +
"dsmcRhoN_",
802 this->
name() +
"linearKE_",
815 this->
name() +
"internalE_",
828 this->
name() +
"iDof_",
841 this->
name() +
"momentum_",
851 rndGen_(Pstream::myProcNo()),
884 binaryCollisionModel_(),
885 wallInteractionModel_(),
886 inflowBoundaryModel_()
890 initialise(dsmcInitialiseDict);
896 template<
class ParcelType>
903 template<
class ParcelType>
906 typename ParcelType::trackingData td(*
this);
913 this->dumpParticlePositions();
917 this->inflowBoundary().inflow();
923 buildCellOccupancy();
933 template<
class ParcelType>
936 label nDSMCParticles = this->size();
939 scalar nMol = nDSMCParticles*nParticle_;
941 vector linearMomentum = linearMomentumOfSystem();
944 scalar linearKineticEnergy = linearKineticEnergyOfSystem();
947 scalar internalEnergy = internalEnergyOfSystem();
951 <<
" Number of dsmc particles = " 957 Info<<
" Number of molecules = " 959 <<
" Mass in system = " 961 <<
" Average linear momentum = " 962 << linearMomentum/nMol <<
nl 963 <<
" |Average linear momentum| = " 964 <<
mag(linearMomentum)/nMol <<
nl 965 <<
" Average linear kinetic energy = " 966 << linearKineticEnergy/nMol <<
nl 967 <<
" Average internal energy = " 968 << internalEnergy/nMol <<
nl 969 <<
" Average total energy = " 970 << (internalEnergy + linearKineticEnergy)/nMol
976 template<
class ParcelType>
985 *rndGen_.GaussNormal<
vector>();
989 template<
class ParcelType>
1005 -
log(rndGen_.sample01<scalar>())
1011 const scalar a = 0.5*iDof - 1;
1012 scalar energyRatio = 0;
1017 energyRatio = 10*rndGen_.sample01<scalar>();
1018 P =
pow((energyRatio/a), a)*
exp(a - energyRatio);
1019 }
while (P < rndGen_.sample01<scalar>());
1025 template<
class ParcelType>
1030 this->db().time().
path()/
"parcelPositions_" 1031 + this->
name() +
"_" 1032 + this->db().time().
timeName() +
".obj" 1035 for (
const ParcelType&
p : *
this)
1037 pObj<<
"v " <<
p.position().x()
1038 <<
' ' <<
p.position().y()
1039 <<
' ' <<
p.position().z()
1047 template<
class ParcelType>
1053 cellOccupancy_.setSize(mesh_.nCells());
1054 buildCellOccupancy();
1057 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.
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 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)