From aaa3d026e8530b4de39a0c1fc58a5fbbeecbc942 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        | 144 +++++++++++++-----
 .../random/cachedRandom/cachedRandom.H        |  53 +++++--
 .../cachedRandom/cachedRandomTemplates.C      |  37 ++++-
 .../StochasticDispersionRAS.C                 |  16 +-
 5 files changed, 179 insertions(+), 75 deletions(-)

diff --git a/src/OpenFOAM/primitives/random/Random/Random.C b/src/OpenFOAM/primitives/random/Random/Random.C
index c0c361e53c..3fa14012c1 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 bbbc0dd344..c7e65ac558 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,13 +121,6 @@ Foam::cachedRandom::~cachedRandom()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-template<>
-Foam::label Foam::cachedRandom::sample01()
-{
-    return round(scalar01());
-}
-
-
 template<>
 Foam::scalar Foam::cachedRandom::sample01()
 {
@@ -127,9 +129,44 @@ Foam::scalar Foam::cachedRandom::sample01()
 
 
 template<>
-Foam::label Foam::cachedRandom::position(const label& start, const label& end)
+Foam::label Foam::cachedRandom::sample01()
 {
-    return start + round(scalar01()*(end - start));
+    return round(scalar01());
+}
+
+
+template<>
+Foam::scalar Foam::cachedRandom::GaussNormal()
+{
+    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::GaussNormal()
+{
+    return round(GaussNormal<scalar>());
 }
 
 
@@ -145,18 +182,9 @@ Foam::scalar Foam::cachedRandom::position
 
 
 template<>
-Foam::label Foam::cachedRandom::globalSample01()
+Foam::label Foam::cachedRandom::position(const label& start, const label& end)
 {
-    scalar value = -GREAT;
-
-    if (Pstream::master())
-    {
-        value = scalar01();
-    }
-
-    reduce(value, maxOp<scalar>());
-
-    return round(value);
+    return start + round(scalar01()*(end - start));
 }
 
 
@@ -170,29 +198,57 @@ Foam::scalar Foam::cachedRandom::globalSample01()
         value = scalar01();
     }
 
-    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::globalSample01()
 {
-    label value = labelMin;
+    scalar value = -GREAT;
 
     if (Pstream::master())
     {
-        value = round(scalar01()*(end - start));
+        value = scalar01();
     }
 
-    reduce(value, maxOp<label>());
+    Pstream::scatter(value);
 
-    return start + value;
+    return round(value);
+}
+
+
+template<>
+Foam::scalar Foam::cachedRandom::globalGaussNormal()
+{
+    scalar value = -GREAT;
+
+    if (Pstream::master())
+    {
+        value = GaussNormal<scalar>();
+    }
+
+    Pstream::scatter(value);
+
+    return value;
+}
+
+
+template<>
+Foam::label Foam::cachedRandom::globalGaussNormal()
+{
+    scalar value = -GREAT;
+
+    if (Pstream::master())
+    {
+        value = GaussNormal<scalar>();
+    }
+
+    Pstream::scatter(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 ec8487fac1..46f23adfe5 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<>
-label cachedRandom::sample01<label>();
-
 template<>
 scalar cachedRandom::sample01<scalar>();
 
 template<>
-label cachedRandom::position<label>(const label& start, const label& end);
+label cachedRandom::sample01<label>();
+
+template<>
+scalar cachedRandom::GaussNormal<scalar>();
+
+template<>
+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 62bb6d4a6a..6684a64cca 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 d2ff74223f..554967aa99 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