ParticleCollector.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2012-2023 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include "ParticleCollector.H"
27 #include "Pstream.H"
28 #include "surfaceWriter.H"
29 #include "unitConversion.H"
30 #include "Random.H"
31 #include "triangle.H"
32 #include "cloud.H"
33 #include "axesRotation.H"
34 #include "OSspecific.H"
35 
36 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
37 
38 template<class CloudType>
40 (
41  const faceList& faces,
42  const Field<point>& points,
43  const Field<scalar>& area
44 )
45 {
46  // Create the output file if not already created
47  if (log_)
48  {
49  if (debug)
50  {
51  Info<< "Creating output file" << endl;
52  }
53 
54  if (Pstream::master())
55  {
56  // Create directory if does not exist
57  mkDir(this->writeTimeDir());
58 
59  // Open new file at start up
60  outputFilePtr_.reset
61  (
62  new OFstream(this->writeTimeDir()/(type() + ".dat"))
63  );
64 
65  outputFilePtr_()
66  << "# Source : " << type() << nl
67  << "# Bins : " << faces.size() << nl
68  << "# Total area : " << sum(area) << nl;
69 
70  outputFilePtr_()
71  << "# Geometry :" << nl
72  << '#'
73  << tab << "Bin"
74  << tab << "(Centre_x Centre_y Centre_z)"
75  << tab << "Area"
76  << nl;
77 
78  forAll(faces, i)
79  {
80  outputFilePtr_()
81  << '#'
82  << tab << i
83  << tab << faces[i].centre(points)
84  << tab << area[i]
85  << nl;
86  }
87 
88  outputFilePtr_()
89  << '#' << nl
90  << "# Output format:" << nl;
91 
92  forAll(faces, i)
93  {
94  word id = Foam::name(i);
95  word binId = "bin_" + id;
96 
97  outputFilePtr_()
98  << '#'
99  << tab << "Time"
100  << tab << binId
101  << tab << "mass[" << id << "]"
102  << tab << "massFlowRate[" << id << "]"
103  << endl;
104  }
105  }
106  }
107 }
108 
109 
110 template<class CloudType>
112 (
113  const List<Field<point>>& polygons
114 )
115 {
116  mode_ = mtPolygon;
117 
118  label nPoints = 0;
119  forAll(polygons, polyI)
120  {
121  label np = polygons[polyI].size();
122  if (np < 3)
123  {
124  FatalIOErrorInFunction(this->coeffDict())
125  << "polygons must consist of at least 3 points"
126  << exit(FatalIOError);
127  }
128 
129  nPoints += np;
130  }
131 
132  label pointOffset = 0;
133  points_.setSize(nPoints);
134  faces_.setSize(polygons.size());
135  area_.setSize(polygons.size());
136  forAll(faces_, facei)
137  {
138  const Field<point>& polyPoints = polygons[facei];
139  face f(identityMap(polyPoints.size()) + pointOffset);
140  UIndirectList<point>(points_, f) = polyPoints;
141  area_[facei] = f.mag(points_);
142  faces_[facei].transfer(f);
143 
144  pointOffset += polyPoints.size();
145  }
146 }
147 
148 
149 template<class CloudType>
151 {
152  mode_ = mtConcentricCircle;
153 
154  vector origin(this->coeffDict().lookup("origin"));
155 
156  this->coeffDict().lookup("radius") >> radius_;
157  nSector_ = this->coeffDict().template lookup<label>("nSector");
158 
159  label nS = nSector_;
160 
161  vector refDir;
162  if (nSector_ > 1)
163  {
164  refDir = this->coeffDict().lookup("refDir");
165  refDir -= normal_[0]*(normal_[0] & refDir);
166  refDir /= mag(refDir);
167  }
168  else
169  {
170  // Set 4 quadrants for single sector cases
171  nS = 4;
172  refDir = normalised(perpendicular(normal_[0]));
173  }
174 
175  scalar dTheta = 5.0;
176  scalar dThetaSector = 360.0/scalar(nS);
177  label intervalPerSector = max(1, ceil(dThetaSector/dTheta));
178  dTheta = dThetaSector/scalar(intervalPerSector);
179 
180  label nPointPerSector = intervalPerSector + 1;
181 
182  label nPointPerRadius = nS*(nPointPerSector - 1);
183  label nPoint = radius_.size()*nPointPerRadius;
184  label nFace = radius_.size()*nS;
185 
186  // Add origin
187  nPoint++;
188 
189  points_.setSize(nPoint);
190  faces_.setSize(nFace);
191  area_.setSize(nFace);
192 
193  coordSys_ = coordinateSystems::cylindrical
194  (
195  "coordSys",
196  origin,
197  normal_[0],
198  refDir,
199  false
200  );
201 
202  List<label> ptIDs(identityMap(nPointPerRadius));
203 
204  points_[0] = origin;
205 
206  // Points
207  forAll(radius_, radI)
208  {
209  label pointOffset = radI*nPointPerRadius + 1;
210 
211  for (label i = 0; i < nPointPerRadius; i++)
212  {
213  label pI = i + pointOffset;
214  point pCyl(radius_[radI], degToRad(i*dTheta), 0.0);
215  points_[pI] = coordSys_.globalPosition(pCyl);
216  }
217  }
218 
219  // Faces
220  DynamicList<label> facePts(2*nPointPerSector);
221  forAll(radius_, radI)
222  {
223  if (radI == 0)
224  {
225  for (label secI = 0; secI < nS; secI++)
226  {
227  facePts.clear();
228 
229  // Append origin point
230  facePts.append(0);
231 
232  for (label ptI = 0; ptI < nPointPerSector; ptI++)
233  {
234  label i = ptI + secI*(nPointPerSector - 1);
235  label id = ptIDs.fcIndex(i - 1) + 1;
236  facePts.append(id);
237  }
238 
239  label facei = secI + radI*nS;
240 
241  faces_[facei] = face(facePts);
242  area_[facei] = faces_[facei].mag(points_);
243  }
244  }
245  else
246  {
247  for (label secI = 0; secI < nS; secI++)
248  {
249  facePts.clear();
250 
251  label offset = (radI - 1)*nPointPerRadius + 1;
252 
253  for (label ptI = 0; ptI < nPointPerSector; ptI++)
254  {
255  label i = ptI + secI*(nPointPerSector - 1);
256  label id = offset + ptIDs.fcIndex(i - 1);
257  facePts.append(id);
258  }
259  for (label ptI = nPointPerSector-1; ptI >= 0; ptI--)
260  {
261  label i = ptI + secI*(nPointPerSector - 1);
262  label id = offset + nPointPerRadius + ptIDs.fcIndex(i - 1);
263  facePts.append(id);
264  }
265 
266  label facei = secI + radI*nS;
267 
268  faces_[facei] = face(facePts);
269  area_[facei] = faces_[facei].mag(points_);
270  }
271  }
272  }
273 }
274 
275 
276 template<class CloudType>
278 (
279  const point& p1,
280  const point& p2
281 ) const
282 {
283  forAll(faces_, facei)
284  {
285  const label facePoint0 = faces_[facei][0];
286 
287  const point& pf = points_[facePoint0];
288 
289  const scalar d1 = normal_[facei] & (p1 - pf);
290  const scalar d2 = normal_[facei] & (p2 - pf);
291 
292  if (sign(d1) == sign(d2))
293  {
294  // Did not cross polygon plane
295  continue;
296  }
297 
298  // Intersection point
299  const point pIntersect = p1 + (d1/(d1 - d2))*(p2 - p1);
300 
301  // Identify if point is within the bounds of the face. Create triangles
302  // between the intersection point and each edge of the face. If all the
303  // triangle normals point in the same direction as the face normal, then
304  // the particle is within the face. Note that testing for pointHits on
305  // the face's decomposed triangles does not work due to ambiguity along
306  // the diagonals.
307  const face& f = faces_[facei];
308  const vector a = f.area(points_);
309  bool inside = true;
310  for (label i = 0; i < f.size(); ++ i)
311  {
312  const label j = f.fcIndex(i);
313  const triPointRef t(pIntersect, points_[f[i]], points_[f[j]]);
314  if ((a & t.area()) < 0)
315  {
316  inside = false;
317  break;
318  }
319  }
320 
321  // Add to the list of hits
322  if (inside)
323  {
324  hitFaceIDs_.append(facei);
325  }
326  }
327 }
328 
329 
330 template<class CloudType>
332 (
333  const point& p1,
334  const point& p2
335 ) const
336 {
337  label secI = -1;
338 
339  const scalar d1 = normal_[0] & (p1 - coordSys_.origin());
340  const scalar d2 = normal_[0] & (p2 - coordSys_.origin());
341 
342  if (sign(d1) == sign(d2))
343  {
344  // Did not cross plane
345  return;
346  }
347 
348  // Intersection point in cylindrical co-ordinate system
349  const point pCyl = coordSys_.localPosition(p1 + (d1/(d1 - d2))*(p2 - p1));
350 
351  scalar r = pCyl[0];
352 
353  if (r < radius_.last())
354  {
355  label radI = 0;
356  while (r > radius_[radI])
357  {
358  radI++;
359  }
360 
361  if (nSector_ == 1)
362  {
363  secI = 4*radI;
364  }
365  else
366  {
367  scalar theta = pCyl[1] + constant::mathematical::pi;
368 
369  secI =
370  nSector_*radI
371  + floor
372  (
373  scalar(nSector_)*theta/constant::mathematical::twoPi
374  );
375  }
376  }
377 
378  if (secI != -1)
379  {
380  hitFaceIDs_.append(secI);
381  }
382 }
383 
384 
385 // * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
386 
387 template<class CloudType>
389 {
390  const fvMesh& mesh = this->owner().mesh();
391  const Time& time = mesh.time();
392  scalar timeNew = time.value();
393  scalar timeElapsed = timeNew - timeOld_;
394 
395  totalTime_ += timeElapsed;
396 
397  const scalar alpha = (totalTime_ - timeElapsed)/totalTime_;
398  const scalar beta = timeElapsed/totalTime_;
399 
400  forAll(faces_, facei)
401  {
402  massFlowRate_[facei] =
403  alpha*massFlowRate_[facei] + beta*mass_[facei]/timeElapsed;
404  massTotal_[facei] += mass_[facei];
405  }
406 
407  const label proci = Pstream::myProcNo();
408 
409  Info<< type() << " output:" << nl;
410 
411  Field<scalar> faceMassTotal(mass_.size(), 0.0);
412  this->getModelProperty("massTotal", faceMassTotal);
413 
414  Field<scalar> faceMassFlowRate(massFlowRate_.size(), 0.0);
415  this->getModelProperty("massFlowRate", faceMassFlowRate);
416 
417 
418  scalar sumTotalMass = 0.0;
419  scalar sumAverageMFR = 0.0;
420  forAll(faces_, facei)
421  {
422  scalarList allProcMass(Pstream::nProcs());
423  allProcMass[proci] = massTotal_[facei];
424  Pstream::gatherList(allProcMass);
425  faceMassTotal[facei] += sum(allProcMass);
426 
427  scalarList allProcMassFlowRate(Pstream::nProcs());
428  allProcMassFlowRate[proci] = massFlowRate_[facei];
429  Pstream::gatherList(allProcMassFlowRate);
430  faceMassFlowRate[facei] += sum(allProcMassFlowRate);
431 
432  sumTotalMass += faceMassTotal[facei];
433  sumAverageMFR += faceMassFlowRate[facei];
434 
435  if (outputFilePtr_.valid())
436  {
437  outputFilePtr_()
438  << time.name()
439  << tab << facei
440  << tab << faceMassTotal[facei]
441  << tab << faceMassFlowRate[facei]
442  << endl;
443  }
444  }
445 
446  Info<< " sum(total mass) = " << sumTotalMass << nl
447  << " sum(average mass flow rate) = " << sumAverageMFR << nl
448  << endl;
449 
450 
451  if (surfaceFormat_ != "none")
452  {
453  if (Pstream::master())
454  {
456  (
457  surfaceWriter::New(surfaceFormat_, this->coeffDict())
458  );
459 
460  writer->write
461  (
462  this->writeTimeDir(),
463  "collector",
464  points_,
465  faces_,
466  false,
467  "massTotal",
468  faceMassTotal,
469  "massFlowRate",
470  faceMassFlowRate
471  );
472  }
473  }
474 
475 
476  if (resetOnWrite_)
477  {
478  Field<scalar> dummy(faceMassTotal.size(), 0.0);
479  this->setModelProperty("massTotal", dummy);
480  this->setModelProperty("massFlowRate", dummy);
481 
482  timeOld_ = timeNew;
483  totalTime_ = 0.0;
484  }
485  else
486  {
487  this->setModelProperty("massTotal", faceMassTotal);
488  this->setModelProperty("massFlowRate", faceMassFlowRate);
489  }
490 
491  forAll(faces_, facei)
492  {
493  mass_[facei] = 0.0;
494  massTotal_[facei] = 0.0;
495  massFlowRate_[facei] = 0.0;
496  }
497 }
498 
499 
500 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
501 
502 template<class CloudType>
504 (
505  const dictionary& dict,
506  CloudType& owner,
507  const word& modelName
508 )
509 :
510  CloudFunctionObject<CloudType>(dict, owner, modelName, typeName),
511  mode_(mtUnknown),
512  parcelType_(this->coeffDict().lookupOrDefault("parcelType", -1)),
513  removeCollected_(this->coeffDict().lookup("removeCollected")),
514  points_(),
515  faces_(),
516  nSector_(0),
517  radius_(),
518  coordSys_("coordSys", vector::zero, axesRotation(sphericalTensor::I)),
519  normal_(),
520  negateParcelsOppositeNormal_
521  (
522  readBool(this->coeffDict().lookup("negateParcelsOppositeNormal"))
523  ),
524  surfaceFormat_(this->coeffDict().lookup("surfaceFormat")),
525  resetOnWrite_(this->coeffDict().lookup("resetOnWrite")),
526  totalTime_(0.0),
527  mass_(),
528  massTotal_(),
529  massFlowRate_(),
530  log_(this->coeffDict().lookup("log")),
531  outputFilePtr_(),
532  timeOld_(owner.mesh().time().value()),
533  hitFaceIDs_()
534 {
535  normal_ /= mag(normal_);
536 
537  word mode(this->coeffDict().lookup("mode"));
538  if (mode == "polygon")
539  {
540  List<Field<point>> polygons(this->coeffDict().lookup("polygons"));
541 
542  initPolygons(polygons);
543 
544  vector n0(this->coeffDict().lookup("normal"));
545  normal_ = vectorField(faces_.size(), n0);
546  }
547  else if (mode == "polygonWithNormal")
548  {
549  List<Tuple2<Field<point>, vector>> polygonAndNormal
550  (
551  this->coeffDict().lookup("polygons")
552  );
553 
554  List<Field<point>> polygons(polygonAndNormal.size());
555  normal_.setSize(polygonAndNormal.size());
556 
557  forAll(polygons, polyI)
558  {
559  polygons[polyI] = polygonAndNormal[polyI].first();
560  normal_[polyI] = polygonAndNormal[polyI].second();
561  normal_[polyI] /= mag(normal_[polyI]) + rootVSmall;
562  }
563 
564  initPolygons(polygons);
565  }
566  else if (mode == "concentricCircle")
567  {
568  vector n0(this->coeffDict().lookup("normal"));
569  normal_ = vectorField(1, n0);
570 
571  initConcentricCircles();
572  }
573  else
574  {
576  << "Unknown mode " << mode << ". Available options are "
577  << "polygon, polygonWithNormal and concentricCircle"
578  << exit(FatalIOError);
579  }
580 
581  mass_.setSize(faces_.size(), 0.0);
582  massTotal_.setSize(faces_.size(), 0.0);
583  massFlowRate_.setSize(faces_.size(), 0.0);
584 
585  makeLogFile(faces_, points_, area_);
586 }
587 
588 
589 template<class CloudType>
591 (
593 )
594 :
596  mode_(pc.mode_),
597  parcelType_(pc.parcelType_),
598  removeCollected_(pc.removeCollected_),
599  points_(pc.points_),
600  faces_(pc.faces_),
601  nSector_(pc.nSector_),
602  radius_(pc.radius_),
603  coordSys_(pc.coordSys_),
604  normal_(pc.normal_),
605  negateParcelsOppositeNormal_(pc.negateParcelsOppositeNormal_),
606  surfaceFormat_(pc.surfaceFormat_),
607  resetOnWrite_(pc.resetOnWrite_),
608  totalTime_(pc.totalTime_),
609  mass_(pc.mass_),
610  massTotal_(pc.massTotal_),
611  massFlowRate_(pc.massFlowRate_),
612  log_(pc.log_),
613  outputFilePtr_(),
614  timeOld_(0.0),
615  hitFaceIDs_()
616 {}
617 
618 
619 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
620 
621 template<class CloudType>
623 {}
624 
625 
626 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
627 
628 template<class CloudType>
630 (
631  parcelType& p,
632  const scalar dt,
633  const point& position0,
634  bool& keepParticle
635 )
636 {
637  if ((parcelType_ != -1) && (parcelType_ != p.typeId()))
638  {
639  return;
640  }
641 
642  hitFaceIDs_.clear();
643 
644  switch (mode_)
645  {
646  case mtPolygon:
647  {
648  collectParcelPolygon
649  (
650  position0,
651  p.position(this->owner().mesh())
652  );
653  break;
654  }
655  case mtConcentricCircle:
656  {
657  collectParcelConcentricCircles
658  (
659  position0,
660  p.position(this->owner().mesh())
661  );
662  break;
663  }
664  default:
665  {}
666  }
667 
668  forAll(hitFaceIDs_, i)
669  {
670  label facei = hitFaceIDs_[i];
671  scalar m = p.nParticle()*p.mass();
672 
673  if (negateParcelsOppositeNormal_)
674  {
675  scalar Unormal = 0;
676  vector Uhat = p.U();
677  switch (mode_)
678  {
679  case mtPolygon:
680  {
681  Unormal = Uhat & normal_[facei];
682  break;
683  }
684  case mtConcentricCircle:
685  {
686  Unormal = Uhat & normal_[0];
687  break;
688  }
689  default:
690  {}
691  }
692 
693  Uhat /= mag(Uhat) + rootVSmall;
694 
695  if (Unormal < 0)
696  {
697  m = -m;
698  }
699  }
700 
701  // Add mass contribution
702  mass_[facei] += m;
703 
704  if (nSector_ == 1)
705  {
706  mass_[facei + 1] += m;
707  mass_[facei + 2] += m;
708  mass_[facei + 3] += m;
709  }
710 
711  if (removeCollected_)
712  {
713  keepParticle = false;
714  }
715  }
716 }
717 
718 
719 // ************************************************************************* //
Functions used by OpenFOAM that are specific to POSIX compliant operating systems and need to be repl...
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
Templated cloud function object base class.
Templated base class for dsmc cloud.
Definition: DSMCCloud.H:79
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
void setSize(const label)
Reset size of List.
Definition: List.C:281
Function object to collect the parcel mass- and mass flow rate over a set of polygons....
virtual void postMove(parcelType &p, const scalar dt, const point &position0, bool &keepParticle)
Post-move hook.
virtual ~ParticleCollector()
Destructor.
ParticleCollector(const dictionary &dict, CloudType &owner, const word &modelName)
Construct from dictionary.
void write()
Write post-processing info.
Templated 3D SphericalTensor derived from VectorSpace adding construction from 1 component,...
Class to control time during OpenFOAM simulations that is also the top-level objectRegistry.
Definition: Time.H:76
T & first()
Return the first element of the list.
Definition: UListI.H:114
An auto-pointer similar to the STL auto_ptr but with automatic casting to a reference to the type and...
Definition: autoPtr.H:51
A coordinate rotation specified using global axis.
Definition: axesRotation.H:67
A list of keyword definitions, which are a keyword followed by any number of values (e....
Definition: dictionary.H:160
const Type & value() const
Return const reference to value.
const word & name() const
Return const reference to name.
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:101
const Time & time() const
Return the top-level database.
Definition: fvMesh.H:402
const dictionary & coeffDict() const
Return const access to the coefficients dictionary.
Definition: subModelBase.C:128
A class for handling words, derived from string.
Definition: word.H:62
A class representing the concept of 0 used to avoid unnecessary manipulations for objects that are kn...
Definition: zero.H:50
#define FatalIOErrorInFunction(ios)
Report an error message using Foam::FatalIOError.
Definition: error.H:318
const pointField & points
label nPoints
volScalarField alpha(IOobject("alpha", runTime.name(), mesh, IOobject::READ_IF_PRESENT, IOobject::AUTO_WRITE), lambda *max(Ua &U, zeroSensitivity))
autoPtr< CompressibleMomentumTransportModel > New(const volScalarField &rho, const volVectorField &U, const surfaceScalarField &phi, const viscosity &viscosity)
const scalar twoPi(2 *pi)
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
bool readBool(Istream &)
Definition: boolIO.C:66
mode_t mode(const fileName &, const bool checkVariants=true, const bool followLink=true)
Return the file mode.
Definition: POSIX.C:461
dimensionedScalar sign(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.
Definition: label.H:59
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
word name(const bool)
Return a word representation of a bool.
Definition: boolIO.C:39
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
static const char tab
Definition: Ostream.H:259
messageStream Info
static const Identity< scalar > I
Definition: Identity.H:93
vector point
Point is a vector.
Definition: point.H:41
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
Vector< Cmpt > perpendicular(const Vector< Cmpt > &v)
Definition: VectorI.H:166
dimensioned< scalar > mag(const dimensioned< Type > &)
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
Form normalised(const VectorSpace< Form, Cmpt, Ncmpts > &vs)
Definition: VectorSpaceI.H:413
Field< vector > vectorField
Specialisation of Field<T> for vector.
IOerror FatalIOError
triangle< point, const point & > triPointRef
Definition: triPointRef.H:44
labelList identityMap(const label len)
Create identity map (map[i] == i) of given length.
Definition: ListOps.C:104
void offset(label &lst, const label o)
List< face > faceList
Definition: faceListFwd.H:43
static const char nl
Definition: Ostream.H:260
fileType type(const fileName &, const bool checkVariants=true, const bool followLink=true)
Return the file type: directory or file.
Definition: POSIX.C:488
scalar degToRad(const scalar deg)
Conversion from degrees to radians.
labelList f(nPoints)
dictionary dict
volScalarField & p
Unit conversion functions.
mkDir(pdfPath)