diff --git a/src/OpenFOAM/Make/files b/src/OpenFOAM/Make/files index 6bd8cbe7f0..8f35d53c96 100644 --- a/src/OpenFOAM/Make/files +++ b/src/OpenFOAM/Make/files @@ -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 diff --git a/src/OpenFOAM/distributions/distribution/distribution.H b/src/OpenFOAM/distributions/distribution/distribution.H index 5a96bc4c0a..99a580256c 100644 --- a/src/OpenFOAM/distributions/distribution/distribution.H +++ b/src/OpenFOAM/distributions/distribution/distribution.H @@ -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 + Type sample() const; + + //- Sample the distribution into a field virtual tmp 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() 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(); \ + } + 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 + Type sample() const + { + Type value; + for (direction i = 0; i < pTraits::nComponents; ++ i) + { + setComponent(value, i) = + static_cast(*this).sample(); + } + return value; + } + + //- Sample the distribution into a field virtual tmp sample(const label n) const { tmp tResult(new scalarField(n)); diff --git a/src/OpenFOAM/distributions/normal/normal.C b/src/OpenFOAM/distributions/normal/normal.C index b794e38c08..bbadcc29f4 100644 --- a/src/OpenFOAM/distributions/normal/normal.C +++ b/src/OpenFOAM/distributions/normal/normal.C @@ -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::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::distributions::normal::phi @@ -85,7 +65,7 @@ Foam::tmp 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& 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_; } diff --git a/src/OpenFOAM/distributions/normal/normal.H b/src/OpenFOAM/distributions/normal/normal.H index 2b47719340..7f9c9b5aee 100644 --- a/src/OpenFOAM/distributions/normal/normal.H +++ b/src/OpenFOAM/distributions/normal/normal.H @@ -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 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 phi diff --git a/src/OpenFOAM/distributions/standardNormal/standardNormal.C b/src/OpenFOAM/distributions/standardNormal/standardNormal.C new file mode 100644 index 0000000000..e0e0bd8734 --- /dev/null +++ b/src/OpenFOAM/distributions/standardNormal/standardNormal.C @@ -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 . + +\*---------------------------------------------------------------------------*/ + +#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::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(rndGen, 0, 0) +{} + + +Foam::distributions::standardNormal::standardNormal(const standardNormal& d) +: + FieldDistribution(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(); + 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::distributions::standardNormal::x(const label n) const +{ + const scalar x0 = approxErfInv(1 - rootSmall); + const scalar x1 = approxErfInv(rootSmall - 1); + + tmp 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::distributions::standardNormal::PDF +( + const scalarField& x +) const +{ + static const scalar sqrt2Pi = sqrt(2*pi); + return exp(- sqr(x)/2)/sqrt2Pi; +} + + +// ************************************************************************* // diff --git a/src/OpenFOAM/distributions/standardNormal/standardNormal.H b/src/OpenFOAM/distributions/standardNormal/standardNormal.H new file mode 100644 index 0000000000..822c6d2411 --- /dev/null +++ b/src/OpenFOAM/distributions/standardNormal/standardNormal.H @@ -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 . + +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 +{ + // 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 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 clone(const label sampleQ) const + { + return autoPtr(new standardNormal(*this)); + } + + + //- Destructor + virtual ~standardNormal(); + + + // Member Functions + + //- Sample the distribution + virtual scalar sample() const; + + //- Sample the distribution + using FieldDistribution::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 x(const label n) const; + + //- Return the distribution probability density function + virtual tmp PDF(const scalarField& x) const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace distributions +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/OpenFOAM/primitives/randomGenerator/randomGenerator.C b/src/OpenFOAM/primitives/randomGenerator/randomGenerator.C index 95cc1f7f7c..d3dbafbc69 100644 --- a/src/OpenFOAM/primitives/randomGenerator/randomGenerator.C +++ b/src/OpenFOAM/primitives/randomGenerator/randomGenerator.C @@ -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; diff --git a/src/OpenFOAM/primitives/randomGenerator/randomGenerator.H b/src/OpenFOAM/primitives/randomGenerator/randomGenerator.H index 83c21b32ed..963b552bbe 100644 --- a/src/OpenFOAM/primitives/randomGenerator/randomGenerator.H +++ b/src/OpenFOAM/primitives/randomGenerator/randomGenerator.H @@ -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 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 - 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(); - // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/OpenFOAM/primitives/randomGenerator/randomGeneratorI.H b/src/OpenFOAM/primitives/randomGenerator/randomGeneratorI.H index 6b953726f7..5141d4844e 100644 --- a/src/OpenFOAM/primitives/randomGenerator/randomGeneratorI.H +++ b/src/OpenFOAM/primitives/randomGenerator/randomGeneratorI.H @@ -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 -inline Type Foam::randomGenerator::sampleNormal() -{ - Type value; - - for (direction i = 0; i < pTraits::nComponents; ++ i) - { - value.component(i) = scalarNormal(); - } - - return value; -} - - -template<> -inline Foam::scalar Foam::randomGenerator::sampleNormal() -{ - return scalarNormal(); -} - - template inline void Foam::randomGenerator::permute(Container& l) { diff --git a/src/lagrangian/DSMC/clouds/Templates/DSMCCloud/DSMCCloud.C b/src/lagrangian/DSMC/clouds/Templates/DSMCCloud/DSMCCloud.C index 3d35beb48a..b2a842321e 100644 --- a/src/lagrangian/DSMC/clouds/Templates/DSMCCloud/DSMCCloud.C +++ b/src/lagrangian/DSMC/clouds/Templates/DSMCCloud/DSMCCloud.C @@ -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::equipartitionLinearVelocity { return sqrt(physicoChemical::k.value()*temperature/mass) - *rndGen_.sampleNormal(); + *distributions::standardNormal(rndGen_).sample(); } diff --git a/src/lagrangian/DSMC/submodels/InflowBoundaryModel/FreeStream/FreeStream.C b/src/lagrangian/DSMC/submodels/InflowBoundaryModel/FreeStream/FreeStream.C index e424abeef9..7a2e4cff85 100644 --- a/src/lagrangian/DSMC/submodels/InflowBoundaryModel/FreeStream/FreeStream.C +++ b/src/lagrangian/DSMC/submodels/InflowBoundaryModel/FreeStream/FreeStream.C @@ -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::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::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; diff --git a/src/lagrangian/DSMC/submodels/WallInteractionModel/MaxwellianThermal/MaxwellianThermal.C b/src/lagrangian/DSMC/submodels/WallInteractionModel/MaxwellianThermal/MaxwellianThermal.C index 4aca443a2b..4faa887396 100644 --- a/src/lagrangian/DSMC/submodels/WallInteractionModel/MaxwellianThermal/MaxwellianThermal.C +++ b/src/lagrangian/DSMC/submodels/WallInteractionModel/MaxwellianThermal/MaxwellianThermal.C @@ -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::correct CloudType& cloud(this->owner()); randomGenerator& rndGen(cloud.rndGen()); + distributions::standardNormal stdNormal(rndGen); while (mag(Ut) < small) { @@ -115,8 +117,8 @@ void Foam::MaxwellianThermal::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 ); diff --git a/src/lagrangian/DSMC/submodels/WallInteractionModel/MixedDiffuseSpecular/MixedDiffuseSpecular.C b/src/lagrangian/DSMC/submodels/WallInteractionModel/MixedDiffuseSpecular/MixedDiffuseSpecular.C index e8b3e55e99..00292f46c5 100644 --- a/src/lagrangian/DSMC/submodels/WallInteractionModel/MixedDiffuseSpecular/MixedDiffuseSpecular.C +++ b/src/lagrangian/DSMC/submodels/WallInteractionModel/MixedDiffuseSpecular/MixedDiffuseSpecular.C @@ -24,6 +24,7 @@ License \*---------------------------------------------------------------------------*/ #include "MixedDiffuseSpecular.H" +#include "standardNormal.H" // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // @@ -79,6 +80,7 @@ void Foam::MixedDiffuseSpecular::correct CloudType& cloud(this->owner()); randomGenerator& rndGen(cloud.rndGen()); + distributions::standardNormal stdNormal(rndGen); if (diffuseFraction_ > rndGen.scalar01()) { @@ -120,8 +122,8 @@ void Foam::MixedDiffuseSpecular::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 ); diff --git a/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudI.H b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudI.H index d5cecb7d78..5725aecdb2 100644 --- a/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudI.H +++ b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudI.H @@ -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(); + *distributions::standardNormal(rndGen_).sample(); } @@ -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() ); } } diff --git a/src/lagrangian/parcel/submodels/Momentum/DispersionModel/GradientDispersionRAS/GradientDispersionRAS.C b/src/lagrangian/parcel/submodels/Momentum/DispersionModel/GradientDispersionRAS/GradientDispersionRAS.C index 3ae06de325..94d7456800 100644 --- a/src/lagrangian/parcel/submodels/Momentum/DispersionModel/GradientDispersionRAS/GradientDispersionRAS.C +++ b/src/lagrangian/parcel/submodels/Momentum/DispersionModel/GradientDispersionRAS/GradientDispersionRAS.C @@ -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::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::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; diff --git a/src/lagrangian/parcel/submodels/Momentum/DispersionModel/StochasticDispersionRAS/StochasticDispersionRAS.C b/src/lagrangian/parcel/submodels/Momentum/DispersionModel/StochasticDispersionRAS/StochasticDispersionRAS.C index 5a96968f32..42e24bd7e4 100644 --- a/src/lagrangian/parcel/submodels/Momentum/DispersionModel/StochasticDispersionRAS/StochasticDispersionRAS.C +++ b/src/lagrangian/parcel/submodels/Momentum/DispersionModel/StochasticDispersionRAS/StochasticDispersionRAS.C @@ -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::update ) { randomGenerator& rndGen = this->owner().rndGen(); + distributions::standardNormal stdNormal(rndGen); const scalar cps = 0.16432; @@ -105,7 +107,7 @@ Foam::vector Foam::StochasticDispersionRAS::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 diff --git a/src/lagrangian/parcel/submodels/Thermodynamic/ParticleForces/BrownianMotion/BrownianMotionForce.C b/src/lagrangian/parcel/submodels/Thermodynamic/ParticleForces/BrownianMotion/BrownianMotionForce.C index b8dc9e443b..2869f45072 100644 --- a/src/lagrangian/parcel/submodels/Thermodynamic/ParticleForces/BrownianMotion/BrownianMotionForce.C +++ b/src/lagrangian/parcel/submodels/Thermodynamic/ParticleForces/BrownianMotion/BrownianMotionForce.C @@ -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 -Foam::scalar Foam::BrownianMotionForce::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 Foam::tmp Foam::BrownianMotionForce::kModel() const @@ -189,28 +171,18 @@ Foam::forceSuSp Foam::BrownianMotionForce::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(); - // 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(); // 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; } diff --git a/src/lagrangian/parcel/submodels/Thermodynamic/ParticleForces/BrownianMotion/BrownianMotionForce.H b/src/lagrangian/parcel/submodels/Thermodynamic/ParticleForces/BrownianMotion/BrownianMotionForce.H index dc59e059b3..0ee5f1562f 100644 --- a/src/lagrangian/parcel/submodels/Thermodynamic/ParticleForces/BrownianMotion/BrownianMotionForce.H +++ b/src/lagrangian/parcel/submodels/Thermodynamic/ParticleForces/BrownianMotion/BrownianMotionForce.H @@ -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 kModel() const; diff --git a/src/randomProcesses/processes/OUprocess/OUprocess.C b/src/randomProcesses/processes/OUprocess/OUprocess.C index e2a8f05597..68870835a3 100644 --- a/src/randomProcesses/processes/OUprocess/OUprocess.C +++ b/src/randomProcesses/processes/OUprocess/OUprocess.C @@ -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()) ); }