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