Stochastic.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) 2013-2018 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 "Stochastic.H"
27 
28 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
29 
30 template<class CloudType>
32 (
33  const dictionary& dict,
34  CloudType& owner
35 )
36 :
37  IsotropyModel<CloudType>(dict, owner, typeName)
38 {}
39 
40 
41 template<class CloudType>
43 (
44  const Stochastic<CloudType>& cm
45 )
46 :
48 {}
49 
50 
51 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
52 
53 template<class CloudType>
55 {}
56 
57 
58 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
59 
60 template<class CloudType>
62 {
63  static bool isCached = true;
64  static scalar xCached;
65 
66  if (isCached)
67  {
68  isCached = false;
69 
70  return xCached;
71  }
72  else
73  {
74  Random& rndGen = this->owner().rndGen();
75 
76  scalar f, m, x, y;
77 
78  do
79  {
80  x = 2.0*rndGen.sample01<scalar>() - 1.0;
81  y = 2.0*rndGen.sample01<scalar>() - 1.0;
82  m = x*x + y*y;
83  } while (m >= 1.0 || m == 0.0);
84 
85  f = sqrt(-2.0*log(m)/m);
86  xCached = x*f;
87  isCached = true;
88 
89  return y*f;
90  }
91 }
92 
93 
94 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
95 
96 template<class CloudType>
98 {
99  const fvMesh& mesh = this->owner().mesh();
100  const scalar deltaT(this->owner().db().time().deltaTValue());
101  Random& rndGen = this->owner().rndGen();
102 
103  const scalar oneBySqrtThree = sqrt(1.0/3.0);
104 
105  const AveragingMethod<scalar>& volumeAverage =
107  (
108  this->owner().name() + ":volumeAverage"
109  );
110  const AveragingMethod<scalar>& radiusAverage =
112  (
113  this->owner().name() + ":radiusAverage"
114  );
115  const AveragingMethod<vector>& uAverage =
117  (
118  this->owner().name() + ":uAverage"
119  );
120  const AveragingMethod<scalar>& uSqrAverage =
122  (
123  this->owner().name() + ":uSqrAverage"
124  );
125  const AveragingMethod<scalar>& frequencyAverage =
127  (
128  this->owner().name() + ":frequencyAverage"
129  );
130  const AveragingMethod<scalar>& massAverage =
132  (
133  this->owner().name() + ":massAverage"
134  );
135 
136  // calculate time scales and pdf exponent
137  autoPtr<AveragingMethod<scalar>> exponentAveragePtr
138  (
140  (
141  IOobject
142  (
143  this->owner().name() + ":exponentAverage",
144  this->owner().db().time().timeName(),
145  mesh
146  ),
147  this->owner().solution().dict(),
148  mesh
149  )
150  );
151  AveragingMethod<scalar>& exponentAverage = exponentAveragePtr();
152  exponentAverage =
153  exp
154  (
155  - deltaT
156  *this->timeScaleModel_->oneByTau
157  (
158  volumeAverage,
159  radiusAverage,
160  uSqrAverage,
161  frequencyAverage
162  )
163  )();
164 
165  // random sampling
166  forAllIter(typename CloudType, this->owner(), iter)
167  {
168  typename CloudType::parcelType& p = iter();
169  const tetIndices tetIs(p.currentTetIndices());
170 
171  const scalar x = exponentAverage.interpolate(p.coordinates(), tetIs);
172 
173  if (x < rndGen.sample01<scalar>())
174  {
175  const vector r(sampleGauss(), sampleGauss(), sampleGauss());
176 
177  const vector u = uAverage.interpolate(p.coordinates(), tetIs);
178  const scalar uRms =
179  sqrt(max(uSqrAverage.interpolate(p.coordinates(), tetIs), 0.0));
180 
181  p.U() = u + r*uRms*oneBySqrtThree;
182  }
183  }
184 
185  // correction velocity averages
186  autoPtr<AveragingMethod<vector>> uTildeAveragePtr
187  (
189  (
190  IOobject
191  (
192  this->owner().name() + ":uTildeAverage",
193  this->owner().db().time().timeName(),
194  mesh
195  ),
196  this->owner().solution().dict(),
197  mesh
198  )
199  );
200  AveragingMethod<vector>& uTildeAverage = uTildeAveragePtr();
201  forAllIter(typename CloudType, this->owner(), iter)
202  {
203  typename CloudType::parcelType& p = iter();
204  const tetIndices tetIs(p.currentTetIndices());
205  uTildeAverage.add(p.coordinates(), tetIs, p.nParticle()*p.mass()*p.U());
206  }
207  uTildeAverage.average(massAverage);
208 
209  autoPtr<AveragingMethod<scalar>> uTildeSqrAveragePtr
210  (
212  (
213  IOobject
214  (
215  this->owner().name() + ":uTildeSqrAverage",
216  this->owner().db().time().timeName(),
217  mesh
218  ),
219  this->owner().solution().dict(),
220  mesh
221  )
222  );
223  AveragingMethod<scalar>& uTildeSqrAverage = uTildeSqrAveragePtr();
224  forAllIter(typename CloudType, this->owner(), iter)
225  {
226  typename CloudType::parcelType& p = iter();
227  const tetIndices tetIs(p.currentTetIndices());
228  const vector uTilde = uTildeAverage.interpolate(p.coordinates(), tetIs);
229  uTildeSqrAverage.add
230  (
231  p.coordinates(),
232  tetIs,
233  p.nParticle()*p.mass()*magSqr(p.U() - uTilde)
234  );
235  }
236  uTildeSqrAverage.average(massAverage);
237 
238  // conservation correction
239  forAllIter(typename CloudType, this->owner(), iter)
240  {
241  typename CloudType::parcelType& p = iter();
242  const tetIndices tetIs(p.currentTetIndices());
243 
244  const vector u = uAverage.interpolate(p.coordinates(), tetIs);
245  const scalar uRms =
246  sqrt(max(uSqrAverage.interpolate(p.coordinates(), tetIs), 0.0));
247 
248  const vector uTilde = uTildeAverage.interpolate(p.coordinates(), tetIs);
249  const scalar uTildeRms =
250  sqrt
251  (
252  max(uTildeSqrAverage.interpolate(p.coordinates(), tetIs), 0.0)
253  );
254 
255  p.U() = u + (p.U() - uTilde)*uRms/max(uTildeRms, small);
256  }
257 }
258 
259 
260 // ************************************************************************* //
dictionary dict
dimensionedScalar log(const dimensionedScalar &ds)
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 > &)
Stochastic return-to-isotropy model.
Definition: Stochastic.H:67
#define forAllIter(Container, container, iter)
Iterate across all elements in the container object of type.
Definition: UList.H:459
dimensionedScalar sqrt(const dimensionedScalar &ds)
virtual void add(const barycentric &coordinates, const tetIndices &tetIs, const Type &value)=0
Member Functions.
const Type & lookupObject(const word &name) const
Lookup and return the object of the given Type.
Random rndGen(label(0))
scalar y
Type sample01()
Advance the state and return a sample of a given type from a.
Definition: RandomI.H:70
dynamicFvMesh & mesh
dimensionedScalar exp(const dimensionedScalar &ds)
word timeName
Definition: getTimeIndex.H:3
Storage and named access for the indices of a tet which is part of the decomposition of a cell...
Definition: tetIndices.H:81
Stochastic(const dictionary &dict, CloudType &owner)
Construct from components.
Definition: Stochastic.C:32
Random number generator.
Definition: Random.H:57
ParcelType parcelType
Type of parcel the cloud was instantiated for.
Definition: DSMCCloud.H:215
dimensioned< scalar > magSqr(const dimensioned< Type > &)
labelList f(nPoints)
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
Base class for collisional return-to-isotropy models.
Definition: MPPICCloud.H:59
virtual ~Stochastic()
Destructor.
Definition: Stochastic.C:54
virtual Type interpolate(const barycentric &coordinates, const tetIndices &tetIs) const =0
Interpolate.
virtual void calculate()
Member Functions.
Definition: Stochastic.C:97
virtual void average()
Calculate the average.
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:78
An auto-pointer similar to the STL auto_ptr but with automatic casting to a reference to the type and...
Definition: PtrList.H:52
Selector class for relaxation factors, solver type and solution.
Definition: solution.H:48
volScalarField & p
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:92
Templated base class for dsmc cloud.
Definition: DSMCCloud.H:69