From 51087ad0df9e8b40c216bb6c35123865dacc234b Mon Sep 17 00:00:00 2001
From: Henry <Henry>
Date: Thu, 9 Apr 2015 15:45:10 +0100
Subject: [PATCH] cachedRandom: Added GaussNormal functions and applied to
StochasticDispersionRAS Resolves bug-report
http://www.openfoam.org/mantisbt/view.php?id=1650 Thanks to Timo Niemi
---
.../primitives/random/Random/Random.C | 4 +-
.../random/cachedRandom/cachedRandom.C | 118 ++++++++++++++----
.../random/cachedRandom/cachedRandom.H | 49 +++++---
.../cachedRandom/cachedRandomTemplates.C | 37 +++++-
.../StochasticDispersionRAS.C | 16 +--
5 files changed, 164 insertions(+), 60 deletions(-)
diff --git a/src/OpenFOAM/primitives/random/Random/Random.C b/src/OpenFOAM/primitives/random/Random/Random.C
index c0c361e53ca..3fa14012c19 100644
--- a/src/OpenFOAM/primitives/random/Random/Random.C
+++ b/src/OpenFOAM/primitives/random/Random/Random.C
@@ -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;
diff --git a/src/OpenFOAM/primitives/random/cachedRandom/cachedRandom.C b/src/OpenFOAM/primitives/random/cachedRandom/cachedRandom.C
index bbbc0dd3445..c7e65ac558d 100644
--- a/src/OpenFOAM/primitives/random/cachedRandom/cachedRandom.C
+++ b/src/OpenFOAM/primitives/random/cachedRandom/cachedRandom.C
@@ -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;
}
diff --git a/src/OpenFOAM/primitives/random/cachedRandom/cachedRandom.H b/src/OpenFOAM/primitives/random/cachedRandom/cachedRandom.H
index ec8487fac1f..46f23adfe5f 100644
--- a/src/OpenFOAM/primitives/random/cachedRandom/cachedRandom.H
+++ b/src/OpenFOAM/primitives/random/cachedRandom/cachedRandom.H
@@ -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
diff --git a/src/OpenFOAM/primitives/random/cachedRandom/cachedRandomTemplates.C b/src/OpenFOAM/primitives/random/cachedRandom/cachedRandomTemplates.C
index 62bb6d4a6ae..6684a64cca8 100644
--- a/src/OpenFOAM/primitives/random/cachedRandom/cachedRandomTemplates.C
+++ b/src/OpenFOAM/primitives/random/cachedRandom/cachedRandomTemplates.C
@@ -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);
}
diff --git a/src/lagrangian/turbulence/submodels/Kinematic/DispersionModel/StochasticDispersionRAS/StochasticDispersionRAS.C b/src/lagrangian/turbulence/submodels/Kinematic/DispersionModel/StochasticDispersionRAS/StochasticDispersionRAS.C
index d2ff74223f8..554967aa99a 100644
--- a/src/lagrangian/turbulence/submodels/Kinematic/DispersionModel/StochasticDispersionRAS/StochasticDispersionRAS.C
+++ b/src/lagrangian/turbulence/submodels/Kinematic/DispersionModel/StochasticDispersionRAS/StochasticDispersionRAS.C
@@ -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
--
GitLab