Go to the documentation of this file.
33 #include "surfaceInterpolate.H"
44 namespace functionObjects
50 extractEulerianParticles,
72 const label nFaces = fz.size();
79 <<
": Number of faceZone faces (" << allFaces
80 <<
") is less than the number of requested locations ("
87 <<
" faces : " << allFaces <<
nl
100 if (!nInjectorLocations_)
105 const faceZone& fz = mesh_.faceZones()[zoneID_];
114 const label nFaces = fz.size();
115 label nLocations = nInjectorLocations_;
120 scalar fraction = scalar(nFaces)/scalar(nGlobalFaces);
121 nLocations = ceil(fraction*nInjectorLocations_);
124 Pout<<
"nFaces:" << nFaces
125 <<
", nGlobalFaces:" << nGlobalFaces
126 <<
", fraction:" << fraction
127 <<
", nLocations:" << nLocations
145 label nCoarseFaces = 0;
148 fineToCoarseAddr_ = ppa.restrictTopBottomAddressing();
149 nCoarseFaces =
max(fineToCoarseAddr_) + 1;
155 <<
" coarse faces" <<
endl;
192 patchIDs_.setSize(fz.size(), -1);
193 patchFaceIDs_.setSize(fz.size(), -1);
195 label nBlockedFaces = 0;
198 const label facei = fz[localFacei];
200 if (mesh_.isInternalFace(facei))
202 if (alphaf[facei] > alphaThreshold_)
204 blockedFaces[localFacei] =
true;
209 label patchi = mesh_.boundaryMesh().whichPatch(facei);
210 label patchFacei = -1;
214 const auto* cpp = isA<coupledPolyPatch>(pp);
218 patchFacei = (cpp->owner() ? pp.
whichFace(facei) : -1);
220 else if (!isA<emptyPolyPatch>(pp))
225 if (patchFacei == -1)
229 else if (alphafp[patchFacei] > alphaThreshold_)
231 blockedFaces[localFacei] =
true;
234 patchIDs_[localFacei] = patchi;
235 patchFaceIDs_[localFacei] = patchFacei;
251 const label particlei = regionToParticleMap_[regioni];
254 if (
p.faceIHit != -1 && nInjectorLocations_)
257 label coarseFacei = fineToCoarseAddr_[
p.faceIHit];
258 p.faceIHit = globalCoarseFaces_.toGlobal(coarseFacei);
265 if ((pDiameter > minDiameter_) && (pDiameter < maxDiameter_))
270 const point position =
p.VC/(
p.V + ROOTVSMALL);
271 const vector U =
p.VU/(
p.V + ROOTVSMALL);
273 if (nInjectorLocations_)
289 cloud_.addParticle(ip);
291 collectedVolume_ +=
p.V;
294 ++nCollectedParticles_;
299 ++nDiscardedParticles_;
300 discardedVolume_ +=
p.V;
307 const label nNewRegions,
316 labelList oldToNewRegion(particles_.size(), -1);
319 forAll(regionFaceIDs, facei)
321 label newRegioni = regionFaceIDs[facei];
322 label oldRegioni = regions0_[facei];
324 if (newRegioni != -1 && oldRegioni != -1)
328 newToNewRegion[newRegioni] =
329 max(newRegioni, oldToNewRegion[oldRegioni]);
330 oldToNewRegion[oldRegioni] = newRegioni;
339 label nParticle = -1;
342 forAll(newToNewRegion, newRegioni0)
344 label newRegioni = newToNewRegion[newRegioni0];
345 if (newRegions.
insert(newRegioni))
351 newRegionToParticleMap.insert(newRegioni0, nParticle);
359 forAll(oldToNewRegion, oldRegioni)
361 label newRegioni = oldToNewRegion[oldRegioni];
362 if (newRegioni == -1)
366 <<
"Collecting particle from oldRegion:" << oldRegioni
369 collectParticle(time, oldRegioni);
374 label newParticlei = newRegionToParticleMap[newRegioni];
375 label oldParticlei = regionToParticleMap_[oldRegioni];
378 <<
"Combining newRegioni: " << newRegioni
379 <<
"(p:" << newParticlei <<
") and "
380 <<
"oldRegioni: " << oldRegioni
381 <<
"(p:" << oldParticlei <<
")"
384 newParticles[newParticlei] =
387 newParticles[newParticlei],
388 particles_[oldParticlei]
394 particles_.transfer(newParticles);
395 regionToParticleMap_ = newRegionToParticleMap;
399 regions0_ = regionFaceIDs;
416 const scalar deltaT = mesh_.time().deltaTValue();
417 const pointField& faceCentres = mesh_.faceCentres();
419 forAll(regionFaceIDs, localFacei)
421 const label newRegioni = regionFaceIDs[localFacei];
422 if (newRegioni != -1)
424 const label particlei = regionToParticleMap_[newRegioni];
425 const label meshFacei = fz[localFacei];
431 p.faceIHit = localFacei;
438 scalar magPhii =
mag(faceValue(
phi, localFacei, meshFacei));
439 vector Ufi = faceValue(
Uf, localFacei, meshFacei);
440 scalar dV = magPhii*deltaT;
442 p.VC += dV*faceCentres[meshFacei];
460 cloud_(mesh_,
"eulerianParticleCloud"),
466 alphaThreshold_(0.1),
470 nInjectorLocations_(0),
472 globalCoarseFaces_(),
475 regionToParticleMap_(),
476 minDiameter_(ROOTVSMALL),
478 nCollectedParticles_(getProperty<label>(
"nCollectedParticles", 0)),
479 collectedVolume_(getProperty<scalar>(
"collectedVolume", 0)),
480 nDiscardedParticles_(getProperty<label>(
"nDiscardedParticles", 0)),
481 discardedVolume_(getProperty<scalar>(
"discardedVolume", 0))
483 if (mesh_.nSolutionD() != 3)
486 <<
name <<
" function object only applicable to 3-D cases"
505 dict.readEntry(
"faceZone", faceZoneName_);
506 dict.readEntry(
"alpha", alphaName_);
508 dict.readIfPresent(
"alphaThreshold", alphaThreshold_);
509 dict.readIfPresent(
"U", UName_);
510 dict.readIfPresent(
"rho", rhoName_);
511 dict.readIfPresent(
"phi", phiName_);
512 dict.readIfPresent(
"nLocations", nInjectorLocations_);
513 dict.readIfPresent(
"minDiameter", minDiameter_);
514 dict.readIfPresent(
"maxDiameter", maxDiameter_);
518 if (nInjectorLocations_)
541 typeName +
":alphaf",
545 const faceZone& fz = mesh_.faceZones()[zoneID_];
553 boolList blockedFaces(fz.size(),
false);
554 setBlockedFaces(alphaf, fz, blockedFaces);
561 const label nRegionsNew = regionFaceIDs.
nRegions();
569 mesh_.time().value(),
575 accumulateParticleInfo(alphaf, tphi(), regionFaceIDs, fz);
577 Log <<
" Collected particles : " << nCollectedParticles_ <<
nl
578 <<
" Collected volume : " << collectedVolume_ <<
nl
579 <<
" Discarded particles : " << nDiscardedParticles_ <<
nl
580 <<
" Discarded volume : " << discardedVolume_ <<
nl
581 <<
" Particles in progress : " << particles_.size() <<
nl
594 setProperty(
"nCollectedParticles", nCollectedParticles_);
595 setProperty(
"collectedVolume", collectedVolume_);
596 setProperty(
"nDiscardedParticles", nDiscardedParticles_);
597 setProperty(
"discardedVolume", discardedVolume_);
int debug
Static debugging option.
extractEulerianParticles(const word &name, const Time &runTime, const dictionary &dict)
Construct from components.
Class to control time during OpenFOAM simulations that is also the top-level objectRegistry.
label nInjectorLocations_
Number of sample locations to generate.
A class for handling words, derived from Foam::string.
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
virtual bool read(const dictionary &)
Read the field min/max data.
A class for managing temporary objects.
Splits a patch into regions based on a mask field. Result is a globally consistent label list of regi...
const dimensionedScalar alpha
Fine-structure constant: default SI units: [].
virtual void initialiseBins()
Initialise the particle collection bins.
bool read(const char *buf, int32_t &val)
Same as readInt32.
autoPtr< surfaceVectorField > Uf
Primitive patch pair agglomerate method.
static bool master(const label communicator=worldComm)
Am I the master process.
Ostream & endl(Ostream &os)
Add newline and flush stream.
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
Specialization of Foam::functionObject for an Foam::fvMesh, providing a reference to the Foam::fvMesh...
#define forAll(list, i)
Loop across all elements in list.
const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0)
virtual bool read(const dictionary &dict)
Read.
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
virtual void collectParticle(const scalar time, const label regioni)
Collect particles that have passed through the faceZone.
const polyBoundaryMesh & boundaryMesh() const
Return boundaryMesh reference.
virtual void checkFaceZone()
Check that the faceZone is valid.
A subset of mesh faces organised as a primitive patch.
messageStream Info
Information stream (stdout output on master, null elsewhere)
#define DebugInFunction
Report an information message using Foam::Info.
A patch is a list of labels that address the faces in the global face list.
void setSize(const label n)
Alias for resize()
virtual void calculateAddressing(const label nRegionsNew, const scalar time, labelList ®ionFaceIDs)
A List with indirect addressing.
const faceZoneMesh & faceZones() const noexcept
Return face zone mesh.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
labelList regions0_
Region indices in faceZone faces from last iteration.
Primarily stores particle properties so that it can be injected at a later time. Note that this store...
virtual bool write()
Write.
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
virtual tmp< surfaceScalarField > phiU() const
Return the volumetric flux.
virtual bool read(const dictionary &dict)
Read optional controls.
Macros for easy insertion into run-time selection tables.
virtual bool write(const token &tok)=0
Write token to stream or otherwise handle it.
const dimensionSet dimMass(1, 0, 0, 0, 0, 0, 0)
static tmp< GeometricField< Type, fvsPatchField, surfaceMesh > > interpolate(const GeometricField< Type, fvPatchField, volMesh > &tvf, const surfaceScalarField &faceFlux, Istream &schemeData)
Interpolate field onto faces using scheme given by Istream.
Calculates a unique integer (label so might not have enough room - 2G max) for processor + local inde...
label findZoneID(const word &zoneName) const
Find zone index by name, return -1 if not found.
wordList names() const
A list of the zone names.
errorManipArg< error, int > exit(error &err, const int errNo=1)
static void listCombineGather(const List< commsStruct > &comms, List< T > &Value, const CombineOp &cop, const int tag, const label comm)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
const word & name() const noexcept
Return the name of this functionObject.
#define DebugInfo
Report an information message using Foam::Info.
label whichFace(const label l) const
Return label of face in patch from global face label.
constexpr scalar pi(M_PI)
label nRegions() const noexcept
Return the global number of regions.
addToRunTimeSelectionTable(functionObject, ObukhovLength, dictionary)
const std::string patch
OpenFOAM patch number as a std::string.
void agglomerate()
Agglomerate patch.
virtual const word & type() const =0
Runtime type information.
static bool & parRun() noexcept
Test if this a parallel run.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
fileName::Type type(const fileName &name, const bool followLink=true)
Return the file type: DIRECTORY or FILE, normally following symbolic links.
const fvMesh & mesh_
Reference to the fvMesh.
labelList identity(const label len, label start=0)
Create identity map of the given length with (map[i] == i)
bool insert(const Key &key)
Insert a new entry, not overwriting existing entries.
defineTypeNameAndDebug(ObukhovLength, 0)
static const word null
An empty word.
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
static const Vector< scalar > zero
Base class for writing single files from the function objects.
word faceZoneName_
Name of faceZone to sample.
virtual void setBlockedFaces(const surfaceScalarField &alphaf, const faceZone &fz, boolList &blockedFaces)
Set the blocked faces, i.e. where alpha > alpha threshold value.
static void listCombineScatter(const List< commsStruct > &comms, List< T > &Value, const int tag, const label comm)
Scatter data. Reverse of combineGather.
dimensionedScalar cbrt(const dimensionedScalar &ds)
label zoneID_
Index of the faceZone.
virtual bool execute()
Execute.
virtual void accumulateParticleInfo(const surfaceScalarField &alphaf, const surfaceScalarField &phi, const labelList ®ionFaceIDs, const faceZone &fz)
Process latest region information.
const Boundary & boundaryField() const
Return const-reference to the boundary field.
A list of faces which address into the list of points.