KinematicCloudI.H
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 2011-2013 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 "fvmSup.H"
27 #include "SortableList.H"
28 
29 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
30 
31 template<class CloudType>
34 {
35  return cloudCopyPtr_();
36 }
37 
38 
39 template<class CloudType>
41 {
42  return mesh_;
43 }
44 
45 
46 template<class CloudType>
47 inline const Foam::IOdictionary&
49 {
50  return particleProperties_;
51 }
52 
53 
54 template<class CloudType>
55 inline const Foam::IOdictionary&
57 {
58  return outputProperties_;
59 }
60 
61 
62 template<class CloudType>
64 {
65  return outputProperties_;
66 }
67 
68 
69 template<class CloudType>
70 inline const Foam::cloudSolution&
72 {
73  return solution_;
74 }
75 
76 
77 template<class CloudType>
79 {
80  return solution_;
81 }
82 
83 
84 template<class CloudType>
85 inline const typename CloudType::particleType::constantProperties&
87 {
88  return constProps_;
89 }
90 
91 
92 template<class CloudType>
93 inline typename CloudType::particleType::constantProperties&
95 {
96  return constProps_;
97 }
98 
99 
100 template<class CloudType>
101 inline const Foam::dictionary&
103 {
104  return subModelProperties_;
105 }
106 
107 
108 template<class CloudType>
110 {
111  return rho_;
112 }
113 
114 
115 template<class CloudType>
117 {
118  return U_;
119 }
120 
121 
122 template<class CloudType>
124 {
125  return mu_;
126 }
127 
128 
129 template<class CloudType>
131 {
132  return g_;
133 }
134 
135 
136 template<class CloudType>
137 inline Foam::scalar Foam::KinematicCloud<CloudType>::pAmbient() const
138 {
139  return pAmbient_;
140 }
141 
142 
143 template<class CloudType>
145 {
146  return pAmbient_;
147 }
148 
149 
150 template<class CloudType>
151 //inline const typename CloudType::parcelType::forceType&
152 inline const typename Foam::KinematicCloud<CloudType>::forceType&
154 {
155  return forces_;
156 }
157 
158 
159 template<class CloudType>
162 {
163  return forces_;
164 }
165 
166 
167 template<class CloudType>
170 {
171  return functions_;
172 }
173 
174 
175 template<class CloudType>
178 {
179  return injectors_;
180 }
181 
182 
183 template<class CloudType>
186 {
187  return injectors_;
188 }
189 
190 
191 template<class CloudType>
194 {
195  return dispersionModel_;
196 }
197 
198 
199 template<class CloudType>
202 {
203  return dispersionModel_();
204 }
205 
206 
207 template<class CloudType>
210 {
211  return patchInteractionModel_;
212 }
213 
214 
215 template<class CloudType>
218 {
219  return patchInteractionModel_();
220 }
221 
222 
223 template<class CloudType>
226 {
227  return stochasticCollisionModel_();
228 }
229 
230 
231 template<class CloudType>
234 {
235  return stochasticCollisionModel_();
236 }
237 
238 
239 template<class CloudType>
242 {
243  return surfaceFilmModel_();
244 }
245 
246 
247 template<class CloudType>
250 {
251  return surfaceFilmModel_();
252 }
253 
254 
255 template<class CloudType>
256 inline const Foam::vectorIntegrationScheme&
258 {
259  return UIntegrator_;
260 }
261 
262 
263 template<class CloudType>
265 {
266  return this->size();
267 }
268 
269 
270 template<class CloudType>
272 {
273  scalar sysMass = 0.0;
274  forAllConstIter(typename KinematicCloud<CloudType>, *this, iter)
275  {
276  const parcelType& p = iter();
277  sysMass += p.nParticle()*p.mass();
278  }
279 
280  return sysMass;
281 }
282 
283 
284 template<class CloudType>
285 inline Foam::vector
287 {
288  vector linearMomentum(vector::zero);
289 
290  forAllConstIter(typename KinematicCloud<CloudType>, *this, iter)
291  {
292  const parcelType& p = iter();
293 
294  linearMomentum += p.nParticle()*p.mass()*p.U();
295  }
296 
297  return linearMomentum;
298 }
299 
300 
301 template<class CloudType>
302 inline Foam::scalar
304 {
305  scalar linearKineticEnergy = 0.0;
306 
307  forAllConstIter(typename KinematicCloud<CloudType>, *this, iter)
308  {
309  const parcelType& p = iter();
310 
311  linearKineticEnergy += p.nParticle()*0.5*p.mass()*(p.U() & p.U());
312  }
313 
314  return linearKineticEnergy;
315 }
316 
317 
318 template<class CloudType>
319 inline Foam::scalar Foam::KinematicCloud<CloudType>::Dij
320 (
321  const label i,
322  const label j
323 ) const
324 {
325  scalar si = 0.0;
326  scalar sj = 0.0;
327  forAllConstIter(typename KinematicCloud<CloudType>, *this, iter)
328  {
329  const parcelType& p = iter();
330  si += p.nParticle()*pow(p.d(), i);
331  sj += p.nParticle()*pow(p.d(), j);
332  }
333 
334  reduce(si, sumOp<scalar>());
335  reduce(sj, sumOp<scalar>());
336  sj = max(sj, VSMALL);
337 
338  return si/sj;
339 }
340 
341 
342 template<class CloudType>
343 inline Foam::scalar Foam::KinematicCloud<CloudType>::Dmax() const
344 {
345  scalar d = -GREAT;
346  forAllConstIter(typename KinematicCloud<CloudType>, *this, iter)
347  {
348  const parcelType& p = iter();
349  d = max(d, p.d());
350  }
351 
352  reduce(d, maxOp<scalar>());
353 
354  return max(0.0, d);
355 }
356 
357 
358 template<class CloudType>
360 (
361  const scalar fraction
362 ) const
363 {
364  if ((fraction < 0) || (fraction > 1))
365  {
367  (
368  "inline Foam::scalar Foam::KinematicCloud<CloudType>::penetration"
369  "("
370  "const scalar"
371  ") const"
372  )
373  << "fraction should be in the range 0 < fraction < 1"
374  << exit(FatalError);
375  }
376 
377  scalar distance = 0.0;
378 
379  const label nParcel = this->size();
380  globalIndex globalParcels(nParcel);
381  const label nParcelSum = globalParcels.size();
382 
383  if (nParcelSum == 0)
384  {
385  return distance;
386  }
387 
388  // lists of parcels mass and distance from initial injection point
389  List<List<scalar> > procMass(Pstream::nProcs());
390  List<List<scalar> > procDist(Pstream::nProcs());
391 
392  List<scalar>& mass = procMass[Pstream::myProcNo()];
393  List<scalar>& dist = procDist[Pstream::myProcNo()];
394 
395  mass.setSize(nParcel);
396  dist.setSize(nParcel);
397 
398  label i = 0;
399  scalar mSum = 0.0;
400  forAllConstIter(typename KinematicCloud<CloudType>, *this, iter)
401  {
402  const parcelType& p = iter();
403  scalar m = p.nParticle()*p.mass();
404  scalar d = mag(p.position() - p.position0());
405  mSum += m;
406 
407  mass[i] = m;
408  dist[i] = d;
409 
410  i++;
411  }
412 
413  // calculate total mass across all processors
414  reduce(mSum, sumOp<scalar>());
415  Pstream::gatherList(procMass);
416  Pstream::gatherList(procDist);
417 
418  if (Pstream::master())
419  {
420  // flatten the mass lists
421  List<scalar> allMass(nParcelSum, 0.0);
422  SortableList<scalar> allDist(nParcelSum, 0.0);
423  for (label procI = 0; procI < Pstream::nProcs(); procI++)
424  {
426  (
427  allMass,
428  globalParcels.localSize(procI),
429  globalParcels.offset(procI)
430  ).assign(procMass[procI]);
431 
432  // flatten the distance list
434  (
435  allDist,
436  globalParcels.localSize(procI),
437  globalParcels.offset(procI)
438  ).assign(procDist[procI]);
439  }
440 
441  // sort allDist distances into ascending order
442  // note: allMass masses are left unsorted
443  allDist.sort();
444 
445  if (nParcelSum > 1)
446  {
447  const scalar mLimit = fraction*mSum;
448  const labelList& indices = allDist.indices();
449 
450  if (mLimit > (mSum - allMass[indices.last()]))
451  {
452  distance = allDist.last();
453  }
454  else
455  {
456  // assuming that 'fraction' is generally closer to 1 than 0,
457  // loop through in reverse distance order
458  const scalar mThreshold = (1.0 - fraction)*mSum;
459  scalar mCurrent = 0.0;
460  label i0 = 0;
461 
462  forAllReverse(indices, i)
463  {
464  label indI = indices[i];
465 
466  mCurrent += allMass[indI];
467 
468  if (mCurrent > mThreshold)
469  {
470  i0 = i;
471  break;
472  }
473  }
474 
475  if (i0 == indices.size() - 1)
476  {
477  distance = allDist.last();
478  }
479  else
480  {
481  // linearly interpolate to determine distance
482  scalar alpha = (mCurrent - mThreshold)/allMass[indices[i0]];
483  distance =
484  allDist[i0] + alpha*(allDist[i0+1] - allDist[i0]);
485  }
486  }
487  }
488  else
489  {
490  distance = allDist.first();
491  }
492  }
493 
494  Pstream::scatter(distance);
495 
496  return distance;
497 }
498 
499 
500 template<class CloudType>
502 {
503  return rndGen_;
504 }
505 
506 
507 template<class CloudType>
510 {
511  if (cellOccupancyPtr_.empty())
512  {
513  buildCellOccupancy();
514  }
515 
516  return cellOccupancyPtr_();
517 }
518 
519 
520 template<class CloudType>
521 inline const Foam::scalarField&
523 {
524  return cellLengthScale_;
525 }
526 
527 
528 template<class CloudType>
531 {
532  return UTrans_();
533 }
534 
535 
536 template<class CloudType>
539 {
540  return UTrans_();
541 }
542 
543 
544 template<class CloudType>
547 {
548  return UCoeff_();
549 }
550 
551 
552 template<class CloudType>
555 {
556  return UCoeff_();
557 }
558 
559 
560 template<class CloudType>
563 {
564  if (debug)
565  {
566  Info<< "UTrans min/max = " << min(UTrans()).value() << ", "
567  << max(UTrans()).value() << nl
568  << "UCoeff min/max = " << min(UCoeff()).value() << ", "
569  << max(UCoeff()).value() << endl;
570  }
571 
572  if (solution_.coupled())
573  {
574  if (solution_.semiImplicit("U"))
575  {
577  Vdt(mesh_.V()*this->db().time().deltaT());
578 
579  return UTrans()/Vdt - fvm::Sp(UCoeff()/Vdt, U) + UCoeff()/Vdt*U;
580  }
581  else
582  {
584  fvVectorMatrix& fvm = tfvm();
585 
586  fvm.source() = -UTrans()/(this->db().time().deltaT());
587 
588  return tfvm;
589  }
590  }
591 
593 }
594 
595 
596 template<class CloudType>
599 {
600  tmp<volScalarField> tvDotSweep
601  (
602  new volScalarField
603  (
604  IOobject
605  (
606  this->name() + ":vDotSweep",
607  this->db().time().timeName(),
608  this->db(),
609  IOobject::NO_READ,
610  IOobject::NO_WRITE,
611  false
612  ),
613  mesh_,
614  dimensionedScalar("zero", dimless/dimTime, 0.0),
615  zeroGradientFvPatchScalarField::typeName
616  )
617  );
618 
619  volScalarField& vDotSweep = tvDotSweep();
620  forAllConstIter(typename KinematicCloud<CloudType>, *this, iter)
621  {
622  const parcelType& p = iter();
623  const label cellI = p.cell();
624 
625  vDotSweep[cellI] += p.nParticle()*p.areaP()*mag(p.U() - U_[cellI]);
626  }
627 
628  vDotSweep.internalField() /= mesh_.V();
629  vDotSweep.correctBoundaryConditions();
630 
631  return tvDotSweep;
632 }
633 
634 
635 template<class CloudType>
638 {
639  tmp<volScalarField> ttheta
640  (
641  new volScalarField
642  (
643  IOobject
644  (
645  this->name() + ":theta",
646  this->db().time().timeName(),
647  this->db(),
648  IOobject::NO_READ,
649  IOobject::NO_WRITE,
650  false
651  ),
652  mesh_,
653  dimensionedScalar("zero", dimless, 0.0),
654  zeroGradientFvPatchScalarField::typeName
655  )
656  );
657 
658  volScalarField& theta = ttheta();
659  forAllConstIter(typename KinematicCloud<CloudType>, *this, iter)
660  {
661  const parcelType& p = iter();
662  const label cellI = p.cell();
663 
664  theta[cellI] += p.nParticle()*p.volume();
665  }
666 
667  theta.internalField() /= mesh_.V();
669 
670  return ttheta;
671 }
672 
673 
674 template<class CloudType>
677 {
678  tmp<volScalarField> talpha
679  (
680  new volScalarField
681  (
682  IOobject
683  (
684  this->name() + ":alpha",
685  this->db().time().timeName(),
686  this->db(),
687  IOobject::NO_READ,
688  IOobject::NO_WRITE,
689  false
690  ),
691  mesh_,
692  dimensionedScalar("zero", dimless, 0.0)
693  )
694  );
695 
696  scalarField& alpha = talpha().internalField();
697  forAllConstIter(typename KinematicCloud<CloudType>, *this, iter)
698  {
699  const parcelType& p = iter();
700  const label cellI = p.cell();
701 
702  alpha[cellI] += p.nParticle()*p.mass();
703  }
704 
705  alpha /= (mesh_.V()*rho_);
706 
707  return talpha;
708 }
709 
710 
711 template<class CloudType>
714 {
715  tmp<volScalarField> trhoEff
716  (
717  new volScalarField
718  (
719  IOobject
720  (
721  this->name() + ":rhoEff",
722  this->db().time().timeName(),
723  this->db(),
724  IOobject::NO_READ,
725  IOobject::NO_WRITE,
726  false
727  ),
728  mesh_,
729  dimensionedScalar("zero", dimDensity, 0.0)
730  )
731  );
732 
733  scalarField& rhoEff = trhoEff().internalField();
734  forAllConstIter(typename KinematicCloud<CloudType>, *this, iter)
735  {
736  const parcelType& p = iter();
737  const label cellI = p.cell();
738 
739  rhoEff[cellI] += p.nParticle()*p.mass();
740  }
741 
742  rhoEff /= mesh_.V();
743 
744  return trhoEff;
745 }
746 
747 
748 // ************************************************************************* //
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
Templated wall surface film model class.
Top level model for Integration schemes.
const vectorIntegrationScheme & UIntegrator() const
Return reference to velocity integration.
scalar Dij(const label i, const label j) const
Mean diameter Dij.
label size() const
Global sum of localSizes.
Definition: globalIndexI.H:66
scalar massInSystem() const
Total mass in system.
const dimensionSet dimForce
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:78
const labelList & indices() const
Return the list of sorted indices. Updated every sort.
Definition: SortableList.H:90
Field with dimensions and associated with geometry type GeoMesh which is used to size the field and a...
label localSize() const
My local size.
Definition: globalIndexI.H:60
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
Calculate the matrix for implicit and explicit sources.
dimensioned< scalar > mag(const dimensioned< Type > &)
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
functionType & functions()
Optional cloud function objects.
const volScalarField & mu() const
Return carrier gas dynamic viscosity.
scalar linearKineticEnergyOfSystem() const
Total linear kinetic energy in the system.
const scalarField & cellLengthScale() const
Return the cell length scale.
const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0)
Definition: dimensionSets.H:51
T & last()
Return the last element of the list.
Definition: UListI.H:131
const parcelType::constantProperties & constProps() const
Return the constant properties.
const IOdictionary & particleProperties() const
Return particle properties dictionary.
Calculates a unique integer (label so might not have enough room - 2G max) for processor + local inde...
Definition: globalIndex.H:63
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
tmp< GeometricField< Type, fvPatchField, volMesh > > Sp(const volScalarField &sp, const GeometricField< Type, fvPatchField, volMesh > &vf)
Definition: fvcSup.C:67
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:76
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
InternalField & internalField()
Return internal field.
const tmp< volScalarField > vDotSweep() const
Volume swept rate of parcels per cell.
const SurfaceFilmModel< KinematicCloud< CloudType > > & surfaceFilm() const
Return const-access to the surface film model.
A special matrix type and solver, designed for finite volume solutions of scalar equations. Face addressing is used to make all matrix assembly and solution loops vectorise.
Definition: fvPatchField.H:68
tmp< fvVectorMatrix > SU(volVectorField &U) const
Return tmp momentum source term.
const tmp< volScalarField > rhoEff() const
Return the particle effective density field.
const InjectionModelList< KinematicCloud< CloudType > > & injectors() const
Return const access to the injection model.
messageStream Info
const dimensionSet dimDensity
const volScalarField & rho() const
Return carrier gas density.
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:137
void sort()
(stable) sort the list (if changed after construction time)
Definition: SortableList.C:95
Templated base class for kinematic cloud.
T & first()
Return the first element of the list.
Definition: UListI.H:117
scalar Dmax() const
Max diameter.
ParcelType parcelType
Type of parcel the cloud was instantiated for.
Definition: DSMCCloud.H:217
fvMatrix< vector > fvVectorMatrix
Definition: fvMatricesFwd.H:45
scalar distance(const vector &p1, const vector &p2)
Definition: curveTools.C:12
Templated patch interaction model class.
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
Definition: pEqn.H:39
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
Definition: HashTable.H:59
const tmp< volScalarField > theta() const
Return the particle volume fraction field.
static const char nl
Definition: Ostream.H:260
void setSize(const label)
Reset size of List.
Definition: List.C:318
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:91
volScalarField & p
Definition: createFields.H:51
A list that is sorted upon construction or when explicitly requested with the sort() method...
Definition: List.H:65
scalar penetration(const scalar fraction) const
Penetration for fraction [0-1] of the current total mass.
label offset(const label procI) const
Start of procI data.
Definition: globalIndexI.H:48
cachedRandom & rndGen()
Return reference to the random object.
scalar pAmbient() const
Return const-access to the ambient pressure.
const dimensionedScalar alpha
Fine-structure constant: default SI units: [].
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
const DispersionModel< KinematicCloud< CloudType > > & dispersion() const
Return const-access to the dispersion model.
IOdictionary is derived from dictionary and IOobject to give the dictionary automatic IO functionalit...
Definition: IOdictionary.H:53
Templated stochastic collision model class.
const KinematicCloud & cloudCopy() const
Return a reference to the cloud copy.
const dimensionedVector & g() const
Gravity.
const PatchInteractionModel< KinematicCloud< CloudType > > & patchInteraction() const
Return const-access to the patch interaction model.
#define FatalErrorIn(functionName)
Report an error message using Foam::FatalError.
Definition: error.H:314
#define forAllReverse(list, i)
Definition: UList.H:424
vector linearMomentumOfSystem() const
Total linear momentum of the system.
error FatalError
const volVectorField & U() const
Return carrier gas velocity.
void correctBoundaryConditions()
Correct boundary field.
Random number generator.
Definition: cachedRandom.H:63
DimensionedField< scalar, volMesh > & UCoeff()
Return coefficient for carrier phase U equation.
List< DynamicList< parcelType * > > & cellOccupancy()
Return the cell occupancy information for each.
A List obtained as a section of another List.
Definition: SubList.H:53
const dimensionSet dimless(0, 0, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:47
Field< Type > & source()
Definition: fvMatrix.H:291
label nParcels() const
Total number of parcels.
const StochasticCollisionModel< KinematicCloud< CloudType > > & stochasticCollision() const
Return const-access to the stochastic collision model.
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
const cloudSolution & solution() const
Return const access to the solution properties.
A class for managing temporary objects.
Definition: PtrList.H:118
Stores all relevant solution info for cloud.
Definition: cloudSolution.H:51
DimensionedField< vector, volMesh > & UTrans()
Return reference to momentum source.
const dictionary & subModelProperties() const
Return reference to the sub-models dictionary.
U
Definition: pEqn.H:82
word timeName
Definition: getTimeIndex.H:3
const tmp< volScalarField > alpha() const
Return the particle mass fraction field.
const forceType & forces() const
Optional particle forces.
const fvMesh & mesh() const
Return reference to the mesh.
List of injection models.
const IOdictionary & outputProperties() const
Return output properties dictionary.