normal.C
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-2015 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 "normal.H"
28 #include "mathematicalConstants.H"
29 
30 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
31 
32 namespace Foam
33 {
34  namespace distributionModels
35  {
36  defineTypeNameAndDebug(normal, 0);
37  addToRunTimeSelectionTable(distributionModel, normal, dictionary);
38  }
39 }
40 
41 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
42 
44 (
45  const dictionary& dict,
46  cachedRandom& rndGen
47 )
48 :
49  distributionModel(typeName, dict, rndGen),
50  minValue_(readScalar(distributionModelDict_.lookup("minValue"))),
51  maxValue_(readScalar(distributionModelDict_.lookup("maxValue"))),
52  expectation_(readScalar(distributionModelDict_.lookup("expectation"))),
53  variance_(readScalar(distributionModelDict_.lookup("variance"))),
54  a_(0.147)
55 {
56  if (minValue_ < 0)
57  {
59  << "Minimum value must be greater than zero. "
60  << "Supplied minValue = " << minValue_
61  << abort(FatalError);
62  }
63 
64  if (maxValue_ < minValue_)
65  {
67  << "Maximum value is smaller than the minimum value:"
68  << " maxValue = " << maxValue_ << ", minValue = " << minValue_
69  << abort(FatalError);
70  }
71 }
72 
73 
75 :
77  minValue_(p.minValue_),
78  maxValue_(p.maxValue_),
79  expectation_(p.expectation_),
80  variance_(p.variance_),
81  a_(p.a_)
82 {}
83 
84 
85 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
86 
88 {}
89 
90 
91 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
92 
94 {
95 
96  scalar a = erf((minValue_ - expectation_)/variance_);
97  scalar b = erf((maxValue_ - expectation_)/variance_);
98 
99  scalar y = rndGen_.sample01<scalar>();
100  scalar x = erfInv(y*(b - a) + a)*variance_ + expectation_;
101 
102  // Note: numerical approximation of the inverse function yields slight
103  // inaccuracies
104 
105  x = min(max(x, minValue_), maxValue_);
106 
107  return x;
108 }
109 
110 
112 {
113  return minValue_;
114 }
115 
116 
118 {
119  return maxValue_;
120 }
121 
122 
124 {
125  return expectation_;
126 }
127 
128 
129 Foam::scalar Foam::distributionModels::normal::erfInv(const scalar y) const
130 {
131  scalar k = 2.0/(constant::mathematical::pi*a_) + 0.5*log(1.0 - y*y);
132  scalar h = log(1.0 - y*y)/a_;
133  scalar x = sqrt(-k + sqrt(k*k - h));
134  if (y < 0.0)
135  {
136  x *= -1.0;
137  }
138  return x;
139 }
140 
141 
142 // ************************************************************************* //
normal(const dictionary &dict, cachedRandom &rndGen)
Construct from components.
Definition: normal.C:44
dimensionedScalar log(const dimensionedScalar &ds)
virtual scalar sample() const
Sample the distributionModel.
Definition: normal.C:93
error FatalError
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:137
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
virtual scalar maxValue() const
Return the maximum value.
Definition: normal.C:117
dimensionedScalar sqrt(const dimensionedScalar &ds)
defineTypeNameAndDebug(distributionModel, 0)
Random number generator.
Definition: cachedRandom.H:63
label k
Boltzmann constant.
Macros for easy insertion into run-time selection tables.
scalar y
Type sample01()
Return a sample whose components lie in the range 0-1.
const dimensionedScalar b
Wien displacement law constant: default SI units: [m.K].
Definition: createFields.H:27
bool readScalar(const char *buf, doubleScalar &s)
Read whole of buf as a scalar. Return true if succesful.
Definition: doubleScalar.H:63
virtual scalar minValue() const
Return the minimum value.
Definition: normal.C:111
virtual ~normal()
Destructor.
Definition: normal.C:87
errorManip< error > abort(error &err)
Definition: errorManip.H:131
virtual scalar erfInv(const scalar y) const
Definition: normal.C:129
dimensionedScalar erf(const dimensionedScalar &ds)
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
const dimensionedScalar h
Planck constant.
Definition: createFields.H:6
cachedRandom & rndGen_
Reference to the random number generator.
addToRunTimeSelectionTable(distributionModel, exponential, dictionary)
volScalarField & p
virtual scalar meanValue() const
Return the mean value.
Definition: normal.C:123
Namespace for OpenFOAM.