Commit 51087ad0 authored by Henry's avatar Henry
Browse files

cachedRandom: Added GaussNormal functions and applied to StochasticDispersionRAS

Resolves bug-report http://www.openfoam.org/mantisbt/view.php?id=1650
Thanks to Timo Niemi
parent c4340790
......@@ -177,10 +177,10 @@ Foam::scalar Foam::Random::GaussNormal()
{
v1 = 2.0*scalar01() - 1.0;
v2 = 2.0*scalar01() - 1.0;
rsq = v1*v1 + v2*v2;
rsq = sqr(v1) + sqr(v2);
} while (rsq >= 1.0 || rsq == 0.0);
fac = sqrt(-2.0 * log(rsq)/rsq);
fac = sqrt(-2.0*log(rsq)/rsq);
gset = v1*fac;
iset = 1;
......
......@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......@@ -57,7 +57,9 @@ Foam::cachedRandom::cachedRandom(const label seed, const label count)
:
seed_(1),
samples_(0),
sampleI_(-1)
sampleI_(-1),
hasGaussSample_(false),
gaussSample_(0)
{
if (seed > 1)
{
......@@ -84,8 +86,15 @@ Foam::cachedRandom::cachedRandom(const cachedRandom& cr, const bool reset)
:
seed_(cr.seed_),
samples_(cr.samples_),
sampleI_(cr.sampleI_)
sampleI_(cr.sampleI_),
hasGaussSample_(cr.hasGaussSample_),
gaussSample_(cr.gaussSample_)
{
if (reset)
{
hasGaussSample_ = false;
gaussSample_ = 0;
}
if (sampleI_ == -1)
{
WarningIn
......@@ -112,6 +121,13 @@ Foam::cachedRandom::~cachedRandom()
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<>
Foam::scalar Foam::cachedRandom::sample01()
{
return scalar01();
}
template<>
Foam::label Foam::cachedRandom::sample01()
{
......@@ -120,16 +136,37 @@ Foam::label Foam::cachedRandom::sample01()
template<>
Foam::scalar Foam::cachedRandom::sample01()
Foam::scalar Foam::cachedRandom::GaussNormal()
{
return scalar01();
if (hasGaussSample_)
{
hasGaussSample_ = false;
return gaussSample_;
}
else
{
scalar rsq, v1, v2;
do
{
v1 = 2*scalar01() - 1;
v2 = 2*scalar01() - 1;
rsq = sqr(v1) + sqr(v2);
} while (rsq >= 1 || rsq == 0);
scalar fac = sqrt(-2*log(rsq)/rsq);
gaussSample_ = v1*fac;
hasGaussSample_ = true;
return v2*fac;
}
}
template<>
Foam::label Foam::cachedRandom::position(const label& start, const label& end)
Foam::label Foam::cachedRandom::GaussNormal()
{
return start + round(scalar01()*(end - start));
return round(GaussNormal<scalar>());
}
......@@ -144,6 +181,29 @@ Foam::scalar Foam::cachedRandom::position
}
template<>
Foam::label Foam::cachedRandom::position(const label& start, const label& end)
{
return start + round(scalar01()*(end - start));
}
template<>
Foam::scalar Foam::cachedRandom::globalSample01()
{
scalar value = -GREAT;
if (Pstream::master())
{
value = scalar01();
}
Pstream::scatter(value);
return value;
}
template<>
Foam::label Foam::cachedRandom::globalSample01()
{
......@@ -154,45 +214,41 @@ Foam::label Foam::cachedRandom::globalSample01()
value = scalar01();
}
reduce(value, maxOp<scalar>());
Pstream::scatter(value);
return round(value);
}
template<>
Foam::scalar Foam::cachedRandom::globalSample01()
Foam::scalar Foam::cachedRandom::globalGaussNormal()
{
scalar value = -GREAT;
if (Pstream::master())
{
value = scalar01();
value = GaussNormal<scalar>();
}
reduce(value, maxOp<scalar>());
Pstream::scatter(value);
return value;
}
template<>
Foam::label Foam::cachedRandom::globalPosition
(
const label& start,
const label& end
)
Foam::label Foam::cachedRandom::globalGaussNormal()
{
label value = labelMin;
scalar value = -GREAT;
if (Pstream::master())
{
value = round(scalar01()*(end - start));
value = GaussNormal<scalar>();
}
reduce(value, maxOp<label>());
Pstream::scatter(value);
return start + value;
return round(value);
}
......@@ -210,17 +266,29 @@ Foam::scalar Foam::cachedRandom::globalPosition
value = scalar01()*(end - start);
}
reduce(value, maxOp<scalar>());
Pstream::scatter(value);
return start + value;
}
void Foam::cachedRandom::operator=(const cachedRandom& cr)
template<>
Foam::label Foam::cachedRandom::globalPosition
(
const label& start,
const label& end
)
{
seed_ = cr.seed_;
samples_ = cr.samples_;
sampleI_ = cr.sampleI_;
label value = labelMin;
if (Pstream::master())
{
value = round(scalar01()*(end - start));
}
Pstream::scatter(value);
return start + value;
}
......
......@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......@@ -63,8 +63,6 @@ class cachedRandom;
class cachedRandom
{
private:
// Private data
//- Initial random number seed
......@@ -76,6 +74,12 @@ private:
//- Current sample marker
label sampleI_;
//- Indicator, which tells if there is a stored gaussian sample
bool hasGaussSample_;
//- Stored sample value
scalar gaussSample_;
// Private Member Functions
......@@ -85,7 +89,6 @@ private:
public:
// Constructors
//- Construct given seed and sample count
......@@ -122,10 +125,16 @@ public:
// Evaluation
// Random numbers
//- Return a sample whose components lie in the range 0-1
template<class Type>
Type sample01();
//- Return a sample whose components are normally distributed
// with zero mean and unity variance N(0, 1)
template<class Type>
Type GaussNormal();
//- Return a sample between start and end
template<class Type>
Type position(const Type& start, const Type& end);
......@@ -141,6 +150,11 @@ public:
template<class Type>
Type globalSample01();
//- Return a sample whose components are normally distributed
// with zero mean and unity variance N(0, 1)
template<class Type>
Type globalGaussNormal();
//- Return a sample between start and end
template<class Type>
Type globalPosition(const Type& start, const Type& end);
......@@ -148,26 +162,23 @@ public:
//- Randomise value in the range 0-1
template<class Type>
void globalRandomise01(Type& value);
// Operators
//- Assignment operator
void operator=(const cachedRandom& cr);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// Template specialisations
template<>
scalar cachedRandom::sample01<scalar>();
template<>
label cachedRandom::sample01<label>();
template<>
scalar cachedRandom::sample01<scalar>();
scalar cachedRandom::GaussNormal<scalar>();
template<>
label cachedRandom::position<label>(const label& start, const label& end);
label cachedRandom::GaussNormal<label>();
template<>
scalar cachedRandom::position<scalar>
......@@ -177,13 +188,19 @@ scalar cachedRandom::position<scalar>
);
template<>
label cachedRandom::globalSample01<label>();
label cachedRandom::position<label>(const label& start, const label& end);
template<>
scalar cachedRandom::globalSample01<scalar>();
template<>
label cachedRandom::globalPosition<label>(const label& start, const label& end);
label cachedRandom::globalSample01<label>();
template<>
scalar cachedRandom::globalGaussNormal<scalar>();
template<>
label cachedRandom::globalGaussNormal<label>();
template<>
scalar cachedRandom::globalPosition<scalar>
......@@ -192,6 +209,10 @@ scalar cachedRandom::globalPosition<scalar>
const scalar& end
);
template<>
label cachedRandom::globalPosition<label>(const label& start, const label& end);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
......
......@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......@@ -41,6 +41,19 @@ Type Foam::cachedRandom::sample01()
}
template<class Type>
Type Foam::cachedRandom::GaussNormal()
{
Type value;
for (direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
{
value.component(cmpt) = GaussNormal<scalar>();
}
return value;
}
template<class Type>
Type Foam::cachedRandom::position(const Type& start, const Type& end)
{
......@@ -72,7 +85,23 @@ Type Foam::cachedRandom::globalSample01()
value = sample01<Type>();
}
reduce(value, maxOp<Type>());
Pstream::scatter(value);
return value;
}
template<class Type>
Type Foam::cachedRandom::globalGaussNormal()
{
Type value = -GREAT*pTraits<Type>::one;
if (Pstream::master())
{
value = GaussNormal<Type>();
}
Pstream::scatter(value);
return value;
}
......@@ -88,7 +117,7 @@ Type Foam::cachedRandom::globalPosition(const Type& start, const Type& end)
value = position<Type>(start, end);
}
reduce(value, maxOp<Type>());
Pstream::scatter(value);
return value;
}
......@@ -104,7 +133,7 @@ void Foam::cachedRandom::globalRandomise01(Type& value)
value = sample01<Type>();
}
reduce(value, maxOp<Type>());
Pstream::scatter(value);
}
......
......@@ -105,21 +105,7 @@ Foam::vector Foam::StochasticDispersionRAS<CloudType>::update
const scalar a = sqrt(1 - sqr(u));
const vector dir(a*cos(theta), a*sin(theta), u);
// Numerical Recipes... Ch. 7. Random Numbers...
scalar x1 = 0;
scalar x2 = 0;
scalar rsq = 10;
while ((rsq > 1) || (rsq == 0))
{
x1 = 2*rnd.sample01<scalar>() - 1;
x2 = 2*rnd.sample01<scalar>() - 1;
rsq = x1*x1 + x2*x2;
}
scalar fac = sqrt(-2*log(rsq)/rsq);
fac *= mag(x1);
UTurb = sigma*fac*dir;
UTurb = sigma*mag(rnd.GaussNormal<scalar>())*dir;
}
}
else
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment