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-2024 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  randomGenerator& 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  randomGenerator& 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().name(),
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(mesh));
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().name(),
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(mesh));
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().name(),
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(mesh));
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(mesh));
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 // ************************************************************************* //
scalar y
#define forAllIter(Container, container, iter)
Iterate across all elements in the container object of type.
Definition: UList.H:459
virtual void add(const barycentric &coordinates, const tetIndices &tetIs, const Type &value)=0
Member Functions.
virtual void average()
Calculate the average.
virtual Type interpolate(const barycentric &coordinates, const tetIndices &tetIs) const =0
Interpolate.
Templated base class for dsmc cloud.
Definition: DSMCCloud.H:80
ParcelType parcelType
Type of parcel the cloud was instantiated for.
Definition: DSMCCloud.H:225
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:99
Base class for collisional return-to-isotropy models.
Definition: IsotropyModel.H:60
Stochastic return-to-isotropy model.
Definition: Stochastic.H:70
virtual ~Stochastic()
Destructor.
Definition: Stochastic.C:54
Stochastic(const dictionary &dict, CloudType &owner)
Construct from components.
Definition: Stochastic.C:32
virtual void calculate()
Member Functions.
Definition: Stochastic.C:97
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 list of keyword definitions, which are a keyword followed by any number of values (e....
Definition: dictionary.H:162
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:99
const polyMesh & mesh() const
Return reference to polyMesh.
Definition: fvMesh.H:441
const Type & lookupObject(const word &name) const
Lookup and return the object of the given Type and name.
Random number generator.
Selector class for relaxation factors, solver type and solution.
Definition: solution.H:51
Storage and named access for the indices of a tet which is part of the decomposition of a cell.
Definition: tetIndices.H:82
dimensionedScalar exp(const dimensionedScalar &ds)
word name(const bool)
Return a word representation of a bool.
Definition: boolIO.C:39
dimensionedScalar log(const dimensionedScalar &ds)
dimensionedScalar sqrt(const dimensionedScalar &ds)
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
dimensioned< scalar > magSqr(const dimensioned< Type > &)
labelList f(nPoints)
dictionary dict
randomGenerator rndGen(653213)
volScalarField & p