zhyang-dev/OpenFOAM-12
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

distributions::standardNormal: Added to replace functions in randomGenerator

Browse Source
This commit is contained in:
Will Bainbridge
2024-04-11 11:44:42 +01:00
parent 1d05b224cb
commit d34d9b93cc
19 changed files with 382 additions and 203 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
src/OpenFOAM/Make/files
View File

@ -719,6 +719,7 @@ distributions/fixedValue/fixedValue.C
distributions/multiNormal/multiNormal.C
distributions/normal/normal.C
distributions/RosinRammler/RosinRammler.C
distributions/standardNormal/standardNormal.C
distributions/tabulatedCumulative/tabulatedCumulative.C
distributions/tabulatedDensity/tabulatedDensity.C
distributions/uniform/uniform.C

54
src/OpenFOAM/distributions/distribution/distribution.H
View File

@ -61,6 +61,7 @@ SourceFiles
#include "dictionary.H"
#include "randomGenerator.H"
#include "fieldTypes.H"
#include "scalarField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -109,10 +110,16 @@ protected:
return sampleQ_ - Q_;
}
//- Sample the distribution into components of a primitive type
#define VIRTUAL_SAMPLE_TYPE(Type, nullArg) \
virtual Type CAT(sample, CAPITALIZE(Type))() const = 0;
FOR_ALL_FIELD_TYPES(VIRTUAL_SAMPLE_TYPE);
#undef VIRTUAL_SAMPLE_TYPE
public:
//-Runtime type information
//- Runtime type information
TypeName("distribution");
@ -181,7 +188,11 @@ public:
//- Sample the distribution
virtual scalar sample() const = 0;
//- Sample the distribution
//- Sample the distribution into components of a primitive type
template<class Type>
Type sample() const;
//- Sample the distribution into a field
virtual tmp<scalarField> sample(const label n) const = 0;
//- Return the minimum value
@ -204,6 +215,16 @@ public:
};
#define DISTRIBUTION_TEMPLATED_SAMPLE_TYPE(Type, nullArg) \
template<> \
inline Type Foam::distribution::sample<Type>() const \
{ \
return CAT(sample, CAPITALIZE(Type))(); \
}
FOR_ALL_FIELD_TYPES(DISTRIBUTION_TEMPLATED_SAMPLE_TYPE);
#undef DISTRIBUTION_TEMPLATED_SAMPLE_TYPE
/*---------------------------------------------------------------------------*\
Class FieldDistribution Declaration
\*---------------------------------------------------------------------------*/
@ -213,6 +234,20 @@ class FieldDistribution
:
public Base
{
protected:
// Protected Member Functions
//- Sample the distribution into components of a primitive type
#define VIRTUAL_SAMPLE_TYPE(Type, nullArg) \
virtual Type CAT(sample, CAPITALIZE(Type))() const \
{ \
return sample<Type>(); \
}
FOR_ALL_FIELD_TYPES(VIRTUAL_SAMPLE_TYPE);
#undef VIRTUAL_SAMPLE_TYPE
public:
// Constructors
@ -226,7 +261,20 @@ public:
//- Sample the distribution
using Base::sample;
//- Sample the distribution
//- Sample the distribution into components of a primitive type
template<class Type>
Type sample() const
{
Type value;
for (direction i = 0; i < pTraits<Type>::nComponents; ++ i)
{
setComponent(value, i) =
static_cast<const Derived&>(*this).sample();
}
return value;
}
//- Sample the distribution into a field
virtual tmp<scalarField> sample(const label n) const
{
tmp<scalarField> tResult(new scalarField(n));

26
src/OpenFOAM/distributions/normal/normal.C
View File

@ -24,6 +24,7 @@ License
\*---------------------------------------------------------------------------*/
#include "normal.H"
#include "standardNormal.H"
#include "mathematicalConstants.H"
#include "addToRunTimeSelectionTable.H"
@ -39,30 +40,9 @@ namespace distributions
}
const Foam::scalar Foam::distributions::normal::a_ = 0.147;
using namespace Foam::constant::mathematical;
// * * * * * * * * * * * Private Static Member Functions * * * * * * * * * * //
Foam::tmp<Foam::scalarField> Foam::distributions::normal::approxErf
(
const scalarField& x
)
{
return sign(x)*sqrt(1 - exp(-sqr(x)*(4/pi + a_*sqr(x))/(1 + a_*sqr(x))));
}
Foam::scalar Foam::distributions::normal::approxErfInv(const scalar y)
{
const scalar l = log(1 - y*y), b = 2/(pi*a_) + l/2;
return sign(y)*sqrt(-b + sqrt(b*b - l/a_));
}
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
Foam::tmp<Foam::scalarField> Foam::distributions::normal::phi
@ -85,7 +65,7 @@ Foam::tmp<Foam::scalarField> Foam::distributions::normal::Phi
if (q == 0)
{
static const scalar sqrt2 = sqrt(scalar(2));
return (1 + approxErf((x - mu_)/(sigma_*sqrt2)))/2;
return (1 + standardNormal::approxErf((x - mu_)/(sigma_*sqrt2)))/2;
}
else
{
@ -99,7 +79,7 @@ Foam::scalar Foam::distributions::normal::sampleForZeroQ(const scalar s) const
static const scalar sqrt2 = sqrt(scalar(2));
const Pair<scalar>& Phi01 = this->Phi01();
const scalar PhiS = (1 - s)*Phi01[0] + s*Phi01[1];
return approxErfInv(2*PhiS - 1)*sigma_*sqrt2 + mu_;
return standardNormal::approxErfInv(2*PhiS - 1)*sigma_*sqrt2 + mu_;
}

13
src/OpenFOAM/distributions/normal/normal.H
View File

@ -29,8 +29,8 @@ Description
minimum and maximum value, rather than from zero to infinity
\f[
PDF(x) = \frac{1}{\sigma \sqrt{2 \pi}} \\
\exp \left( \frac{1}{2} \left( \frac{x - \mu}{\sigma} \right)^2 \right)
PDF(x) = \frac{1}{\sigma \sqrt{2 \pi}} \exp \\
\left( - \frac{1}{2} \left( \frac{x - \mu}{\sigma} \right)^2 \right)
\f]
Usage
@ -88,18 +88,9 @@ class normal
//- Standard deviation
const scalar sigma_;
//- Constant for approximate error function and inverse error function
static const scalar a_;
// Private Member Functions
//- Approximate error function
static tmp<scalarField> approxErf(const scalarField& x);
//- Approximate error function inverse
static scalar approxErfInv(const scalar y);
//- Return values of the un-normalised PDF for the given size exponent
// and x-coordinates.
virtual tmp<scalarField> phi

142
src/OpenFOAM/distributions/standardNormal/standardNormal.C Normal file
View File

@ -0,0 +1,142 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2024 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "standardNormal.H"
#include "mathematicalConstants.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
namespace distributions
{
defineTypeNameAndDebug(standardNormal, 0);
}
}
const Foam::scalar Foam::distributions::standardNormal::a_ = 0.147;
using namespace Foam::constant::mathematical;
// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
Foam::tmp<Foam::scalarField> Foam::distributions::standardNormal::approxErf
(
const scalarField& x
)
{
return sign(x)*sqrt(1 - exp(-sqr(x)*(4/pi + a_*sqr(x))/(1 + a_*sqr(x))));
}
Foam::scalar Foam::distributions::standardNormal::approxErfInv(const scalar y)
{
const scalar l = log(1 - y*y), b = 2/(pi*a_) + l/2;
return sign(y)*sqrt(-b + sqrt(b*b - l/a_));
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::distributions::standardNormal::standardNormal(randomGenerator& rndGen)
:
FieldDistribution<distribution, standardNormal>(rndGen, 0, 0)
{}
Foam::distributions::standardNormal::standardNormal(const standardNormal& d)
:
FieldDistribution<distribution, standardNormal>(d.rndGen_, 0, 0)
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::distributions::standardNormal::~standardNormal()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::scalar Foam::distributions::standardNormal::sample() const
{
static const scalar sqrt2 = sqrt(scalar(2));
const scalar s = rndGen_.sample01<scalar>();
return approxErfInv(2*s - 1)*sqrt2;
}
Foam::scalar Foam::distributions::standardNormal::min() const
{
return -vGreat;
}
Foam::scalar Foam::distributions::standardNormal::max() const
{
return vGreat;
}
Foam::scalar Foam::distributions::standardNormal::mean() const
{
return 0;
}
Foam::tmp<Foam::scalarField>
Foam::distributions::standardNormal::x(const label n) const
{
const scalar x0 = approxErfInv(1 - rootSmall);
const scalar x1 = approxErfInv(rootSmall - 1);
tmp<scalarField> tResult(new scalarField(n));
scalarField& result = tResult.ref();
forAll(result, i)
{
const scalar f = scalar(i)/(n - 1);
result[i] = (1 - f)*x0 + f*x1;
}
return tResult;
}
Foam::tmp<Foam::scalarField> Foam::distributions::standardNormal::PDF
(
const scalarField& x
) const
{
static const scalar sqrt2Pi = sqrt(2*pi);
return exp(- sqr(x)/2)/sqrt2Pi;
}
// ************************************************************************* //

140
src/OpenFOAM/distributions/standardNormal/standardNormal.H Normal file
View File

@ -0,0 +1,140 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2024 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::distributions::standardNormal
Description
Standard normal distribution. Not selectable.
\f[
PDF(x) = \frac{1}{\sqrt{2 \pi}} \exp \left( - \frac{1}{2} x^2 \right)
\f]
SourceFiles
standardNormal.C
See also
Foam::distribution
\*---------------------------------------------------------------------------*/
#ifndef standardNormal_H
#define standardNormal_H
#include "unintegrable.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace distributions
{
/*---------------------------------------------------------------------------*\
Class standardNormal Declaration
\*---------------------------------------------------------------------------*/
class standardNormal
:
public FieldDistribution<distribution, standardNormal>
{
// Private Data
//- Constant for approximate error function and inverse error function
static const scalar a_;
public:
//- Runtime type information
TypeName("standardNormal");
//- Permit the multiNormal distribution to use private parts of this class
friend class multiNormal;
// Static Member Functions
//- Approximate error function
static tmp<scalarField> approxErf(const scalarField& x);
//- Approximate error function inverse
static scalar approxErfInv(const scalar y);
// Constructors
//- Construct from a random generator
standardNormal(randomGenerator& rndGen);
//- Construct copy
standardNormal(const standardNormal& d);
//- Construct and return a clone
virtual autoPtr<distribution> clone(const label sampleQ) const
{
return autoPtr<distribution>(new standardNormal(*this));
}
//- Destructor
virtual ~standardNormal();
// Member Functions
//- Sample the distribution
virtual scalar sample() const;
//- Sample the distribution
using FieldDistribution<distribution, standardNormal>::sample;
//- Return the minimum value
virtual scalar min() const;
//- Return the maximum value
virtual scalar max() const;
//- Return the mean value
virtual scalar mean() const;
//- Return coordinates to plot across the range of the distribution
virtual tmp<scalarField> x(const label n) const;
//- Return the distribution probability density function
virtual tmp<scalarField> PDF(const scalarField& x) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace distributions
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

63
src/OpenFOAM/primitives/randomGenerator/randomGenerator.C
View File

@ -28,69 +28,6 @@ License
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
Foam::scalar Foam::randomGenerator::scalarNormal()
{
// Proper inversion of the distribution. Slow. Exactly maintains
// the random behaviour of the generator.
/*
using namespace constant::mathematical;
static const scalar sqrtTwo = sqrt(scalar(2));
static const scalar sqrtPiByTwo = sqrt(pi)/2;
static const scalar a = 8*(pi - 3)/(3*pi*(4 - pi));
const scalar x = 2*scalar01() - 1;
const scalar xPos = mag(x);
// Initial approximation
const scalar l = log(1 - sqr(xPos));
const scalar ll = 2/(pi*a) + l/2;
scalar y = sqrt(sqrt(sqr(ll) - l/a) - ll);
// Newton improvement
label n = 0;
while (n < 2)
{
const scalar dt = (erf(y) - xPos)/exp(- y*y)*sqrtPiByTwo;
y -= dt;
n += mag(dt) < rootSmall;
}
return sign(x)*sqrtTwo*y;
*/
// Box-Muller transform. Fast. Uses rejection and caching so the
// random sequence is not guaranteed.
if (scalarNormalStored_)
{
scalarNormalStored_ = false;
return scalarNormalValue_;
}
else
{
scalar x1, x2, rr;
do
{
x1 = 2*scalar01() - 1;
x2 = 2*scalar01() - 1;
rr = sqr(x1) + sqr(x2);
}
while (rr >= 1 || rr == 0);
const scalar f = sqrt(- 2*log(rr)/rr);
scalarNormalValue_ = x1*f;
scalarNormalStored_ = true;
return x2*f;
}
}
Foam::scalar Foam::randomGenerator::globalScalar01()
{
scalar value = - vGreat;

18
src/OpenFOAM/primitives/randomGenerator/randomGenerator.H
View File

@ -75,12 +75,6 @@ class randomGenerator
//- The stored integer
type x_;
//- Is a normal scalar sample stored?
bool scalarNormalStored_;
//- A stored normal scalar sample
scalar scalarNormalValue_;
// Private Member Functions
@ -112,10 +106,6 @@ public:
// distribution between two limits
inline scalar scalarAB(const scalar a, const scalar b);
//- Advance the state and return a scalar sample from a normal
// distribution with mean zero and standard deviation one
scalar scalarNormal();
// Other types
@ -129,11 +119,6 @@ public:
template<class Type>
inline Type sampleAB(const Type& a, const Type& b);
//- Advance the state and return a sample of a given type from a
// normal distribution with mean zero and standard deviation one
template<class Type>
inline Type sampleNormal();
// Global scalars
@ -162,9 +147,6 @@ inline scalar randomGenerator::sampleAB(const scalar& a, const scalar& b);
template<>
inline label randomGenerator::sampleAB(const label& a, const label& b);
template<>
inline scalar randomGenerator::sampleNormal();
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

26
src/OpenFOAM/primitives/randomGenerator/randomGeneratorI.H
View File

@ -30,7 +30,6 @@ License
inline Foam::randomGenerator::type Foam::randomGenerator::sample()
{
x_ = (A*x_ + C) % M;
return x_ >> 17;
}
@ -39,9 +38,7 @@ inline Foam::randomGenerator::type Foam::randomGenerator::sample()
inline Foam::randomGenerator::randomGenerator(const label s)
:
x_((type(s) << 16) + 0x330E),
scalarNormalStored_(false),
scalarNormalValue_(0)
x_((type(s) << 16) + 0x330E)
{}
@ -126,27 +123,6 @@ inline Foam::label Foam::randomGenerator::sampleAB
}
template<class Type>
inline Type Foam::randomGenerator::sampleNormal()
{
Type value;
for (direction i = 0; i < pTraits<Type>::nComponents; ++ i)
{
value.component(i) = scalarNormal();
}
return value;
}
template<>
inline Foam::scalar Foam::randomGenerator::sampleNormal()
{
return scalarNormal();
}
template<class Container>
inline void Foam::randomGenerator::permute(Container& l)
{

5
src/lagrangian/DSMC/clouds/Templates/DSMCCloud/DSMCCloud.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2023 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2024 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -30,6 +30,7 @@ License
#include "constants.H"
#include "zeroGradientFvPatchFields.H"
#include "polyMeshTetDecomposition.H"
#include "standardNormal.H"
using namespace Foam::constant;
@ -1042,7 +1043,7 @@ Foam::vector Foam::DSMCCloud<ParcelType>::equipartitionLinearVelocity
{
return
sqrt(physicoChemical::k.value()*temperature/mass)
*rndGen_.sampleNormal<vector>();
*distributions::standardNormal(rndGen_).sample<vector>();
}

7
src/lagrangian/DSMC/submodels/InflowBoundaryModel/FreeStream/FreeStream.C
View File

@ -27,6 +27,7 @@ License
#include "constants.H"
#include "triPointRef.H"
#include "tetIndices.H"
#include "standardNormal.H"
using namespace Foam::constant::mathematical;
@ -146,6 +147,7 @@ void Foam::FreeStream<CloudType>::inflow()
const scalar deltaT = mesh.time().deltaTValue();
randomGenerator& rndGen(cloud.rndGen());
distributions::standardNormal stdNormal(rndGen);
scalar sqrtPi = sqrt(pi);
@ -388,10 +390,7 @@ void Foam::FreeStream<CloudType>::inflow()
vector U =
sqrt(physicoChemical::k.value()*faceTemperature/mass)
*(
rndGen.scalarNormal()*t1
+ rndGen.scalarNormal()*t2
)
*(stdNormal.sample()*t1 + stdNormal.sample()*t2)
+ (t1 & faceVelocity)*t1
+ (t2 & faceVelocity)*t2
+ mostProbableSpeed*uNormal*n;

6
src/lagrangian/DSMC/submodels/WallInteractionModel/MaxwellianThermal/MaxwellianThermal.C
View File

@ -25,6 +25,7 @@ License
#include "MaxwellianThermal.H"
#include "constants.H"
#include "standardNormal.H"
using namespace Foam::constant;
@ -81,6 +82,7 @@ void Foam::MaxwellianThermal<CloudType>::correct
CloudType& cloud(this->owner());
randomGenerator& rndGen(cloud.rndGen());
distributions::standardNormal stdNormal(rndGen);
while (mag(Ut) < small)
{
@ -115,8 +117,8 @@ void Foam::MaxwellianThermal<CloudType>::correct
U =
sqrt(physicoChemical::k.value()*T/mass)
*(
rndGen.scalarNormal()*tw1
+ rndGen.scalarNormal()*tw2
stdNormal.sample()*tw1
+ stdNormal.sample()*tw2
- sqrt(-2.0*log(max(1 - rndGen.scalar01(), vSmall)))*nw
);

6
src/lagrangian/DSMC/submodels/WallInteractionModel/MixedDiffuseSpecular/MixedDiffuseSpecular.C
View File

@ -24,6 +24,7 @@ License
\*---------------------------------------------------------------------------*/
#include "MixedDiffuseSpecular.H"
#include "standardNormal.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
@ -79,6 +80,7 @@ void Foam::MixedDiffuseSpecular<CloudType>::correct
CloudType& cloud(this->owner());
randomGenerator& rndGen(cloud.rndGen());
distributions::standardNormal stdNormal(rndGen);
if (diffuseFraction_ > rndGen.scalar01())
{
@ -120,8 +122,8 @@ void Foam::MixedDiffuseSpecular<CloudType>::correct
U =
sqrt(physicoChemical::k.value()*T/mass)
*(
rndGen.scalarNormal()*tw1
+ rndGen.scalarNormal()*tw2
stdNormal.sample()*tw1
+ stdNormal.sample()*tw2
- sqrt(-2.0*log(max(1 - rndGen.scalar01(), vSmall)))*nw
);

15
src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudI.H
View File

@ -24,6 +24,7 @@ License
\*---------------------------------------------------------------------------*/
#include "constants.H"
#include "standardNormal.H"
using namespace Foam::constant;
@ -303,7 +304,7 @@ inline Foam::vector Foam::moleculeCloud::equipartitionLinearVelocity
{
return
sqrt(physicoChemical::k.value()*temperature/mass)
*rndGen_.sampleNormal<vector>();
*distributions::standardNormal(rndGen_).sample<vector>();
}
@ -315,22 +316,24 @@ inline Foam::vector Foam::moleculeCloud::equipartitionAngularMomentum
{
scalar sqrtKbT = sqrt(physicoChemical::k.value()*temperature);
distributions::standardNormal stdNormal(rndGen_);
if (cP.linearMolecule())
{
return sqrtKbT*vector
(
0.0,
sqrt(cP.momentOfInertia().yy())*rndGen_.scalarNormal(),
sqrt(cP.momentOfInertia().zz())*rndGen_.scalarNormal()
sqrt(cP.momentOfInertia().yy())*stdNormal.sample(),
sqrt(cP.momentOfInertia().zz())*stdNormal.sample()
);
}
else
{
return sqrtKbT*vector
(
sqrt(cP.momentOfInertia().xx())*rndGen_.scalarNormal(),
sqrt(cP.momentOfInertia().yy())*rndGen_.scalarNormal(),
sqrt(cP.momentOfInertia().zz())*rndGen_.scalarNormal()
sqrt(cP.momentOfInertia().xx())*stdNormal.sample(),
sqrt(cP.momentOfInertia().yy())*stdNormal.sample(),
sqrt(cP.momentOfInertia().zz())*stdNormal.sample()
);
}
}

7
src/lagrangian/parcel/submodels/Momentum/DispersionModel/GradientDispersionRAS/GradientDispersionRAS.C
View File

@ -26,6 +26,7 @@ License
#include "GradientDispersionRAS.H"
#include "demandDrivenData.H"
#include "fvcGrad.H"
#include "standardNormal.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
@ -100,7 +101,7 @@ Foam::vector Foam::GradientDispersionRAS<CloudType>::update
scalar& tTurb
)
{
randomGenerator& rndGen = this->owner().rndGen();
distributions::standardNormal stdNormal(this->owner().rndGen());
const scalar cps = 0.16432;
@ -135,11 +136,11 @@ Foam::vector Foam::GradientDispersionRAS<CloudType>::update
// prevent this we let fac be both negative/positive
if (this->owner().mesh().nSolutionD() == 2)
{
fac = rndGen.scalarNormal();
fac = stdNormal.sample();
}
else
{
fac = mag(rndGen.scalarNormal());
fac = mag(stdNormal.sample());
}
UTurb = sigma*fac*dir;

4
src/lagrangian/parcel/submodels/Momentum/DispersionModel/StochasticDispersionRAS/StochasticDispersionRAS.C
View File

@ -25,6 +25,7 @@ License
#include "StochasticDispersionRAS.H"
#include "constants.H"
#include "standardNormal.H"
using namespace Foam::constant::mathematical;
@ -72,6 +73,7 @@ Foam::vector Foam::StochasticDispersionRAS<CloudType>::update
)
{
randomGenerator& rndGen = this->owner().rndGen();
distributions::standardNormal stdNormal(rndGen);
const scalar cps = 0.16432;
@ -105,7 +107,7 @@ Foam::vector Foam::StochasticDispersionRAS<CloudType>::update
const scalar a = sqrt(1 - sqr(u));
const vector dir(a*cos(theta), a*sin(theta), u);
UTurb = sigma*mag(rndGen.scalarNormal())*dir;
UTurb = sigma*mag(stdNormal.sample())*dir;
}
}
else

40
src/lagrangian/parcel/submodels/Thermodynamic/ParticleForces/BrownianMotion/BrownianMotionForce.C
View File

@ -28,30 +28,12 @@ License
#include "fundamentalConstants.H"
#include "demandDrivenData.H"
#include "momentumTransportModel.H"
#include "standardNormal.H"
using namespace Foam::constant;
// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
template<class CloudType>
Foam::scalar Foam::BrownianMotionForce<CloudType>::erfInv(const scalar y) const
{
const scalar a = 0.147;
scalar k = 2.0/(mathematical::pi*a) + 0.5*log(1.0 - y*y);
scalar h = log(1.0 - y*y)/a;
scalar x = sqrt(-k + sqrt(k*k - h));
if (y < 0.0)
{
return -x;
}
else
{
return x;
}
}
template<class CloudType>
Foam::tmp<Foam::volScalarField>
Foam::BrownianMotionForce<CloudType>::kModel() const
@ -189,28 +171,18 @@ Foam::forceSuSp Foam::BrownianMotionForce<CloudType>::calcCoupled
f = mass*sqrt(mathematical::pi*s0/dt);
}
randomGenerator& rndGen = this->owner().rndGen();
distributions::standardNormal stdNormal(rndGen);
// To generate a cubic distribution (3 independent directions) :
// const scalar sqrt2 = sqrt(2.0);
// for (direction dir = 0; dir < vector::nComponents; dir++)
// {
// const scalar x = rndGen_.sample01<scalar>();
// const scalar eta = sqrt2*erfInv(2*x - 1.0);
// value.Su()[dir] = f*eta;
// }
// To generate a cubic distribution (i.e., 3 independent directions):
// value.Su() = f*stdNormal.sample<vector>();
// To generate a spherical distribution:
randomGenerator& rndGen = this->owner().rndGen();
const scalar theta = rndGen.scalar01()*twoPi;
const scalar u = 2*rndGen.scalar01() - 1;
const scalar a = sqrt(1 - sqr(u));
const vector dir(a*cos(theta), a*sin(theta), u);
value.Su() = f*mag(rndGen.scalarNormal())*dir;
value.Su() = f*mag(stdNormal.sample())*dir;
return value;
}

3
src/lagrangian/parcel/submodels/Thermodynamic/ParticleForces/BrownianMotion/BrownianMotionForce.H
View File

@ -82,9 +82,6 @@ class BrownianMotionForce
// Private Member Functions
//- Inverse erf for Gaussian distribution
scalar erfInv(const scalar y) const;
//- Return the k field from the turbulence model
tmp<volScalarField> kModel() const;

9
src/randomProcesses/processes/OUprocess/OUprocess.C
View File

@ -25,6 +25,7 @@ License
#include "OUprocess.H"
#include "Kmesh.H"
#include "standardNormal.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -35,11 +36,13 @@ namespace Foam
complexVector OUprocess::WeinerProcess(const scalar deltaT) const
{
distributions::standardNormal stdNormal(rndGen_);
return sqrt(deltaT)*complexVector
(
complex(rndGen_.scalarNormal(), rndGen_.scalarNormal()),
complex(rndGen_.scalarNormal(), rndGen_.scalarNormal()),
complex(rndGen_.scalarNormal(), rndGen_.scalarNormal())
complex(stdNormal.sample(), stdNormal.sample()),
complex(stdNormal.sample(), stdNormal.sample()),
complex(stdNormal.sample(), stdNormal.sample())
);
}