Merge branch 'feature-distributionModels' into 'develop'
Improve the Lagrangian distribution models See merge request Development/openfoam!426
This commit is contained in:
Andrew Heather
commit
c307c4abd2
@ -6,6 +6,7 @@
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
Copyright (C) 2011-2020 OpenFOAM Foundation
|
||||
Copyright (C) 2021 OpenCFD Ltd.
|
||||
-------------------------------------------------------------------------------
|
||||
License
|
||||
This file is part of OpenFOAM.
|
||||
@ -26,6 +27,7 @@ License
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
#include "RosinRammler.H"
|
||||
#include "MathFunctions.H"
|
||||
#include "addToRunTimeSelectionTable.H"
|
||||
|
||||
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
|
||||
@ -48,11 +50,29 @@ Foam::distributionModels::RosinRammler::RosinRammler
|
||||
)
|
||||
:
|
||||
distributionModel(typeName, dict, rndGen),
|
||||
minValue_(distributionModelDict_.get<scalar>("minValue")),
|
||||
maxValue_(distributionModelDict_.get<scalar>("maxValue")),
|
||||
d_(distributionModelDict_.get<scalar>("d")),
|
||||
lambda_(distributionModelDict_.getCompat<scalar>("lambda", {{"d", 2112}})),
|
||||
n_(distributionModelDict_.get<scalar>("n"))
|
||||
{
|
||||
const word parcelBasisType =
|
||||
dict.getOrDefault<word>("parcelBasisType", "none");
|
||||
|
||||
if (parcelBasisType == "mass")
|
||||
{
|
||||
WarningInFunction
|
||||
<< "Selected parcel basis type: " << parcelBasisType << nl
|
||||
<< " Please consider to use massRosinRammler distribution model"
|
||||
<< endl;
|
||||
}
|
||||
|
||||
if (lambda_ < VSMALL || n_ < VSMALL)
|
||||
{
|
||||
FatalErrorInFunction
|
||||
<< "Scale/Shape parameter cannot be equal to or less than zero:"
|
||||
<< " lambda = " << lambda_
|
||||
<< " n = " << n_
|
||||
<< exit(FatalError);
|
||||
}
|
||||
|
||||
check();
|
||||
}
|
||||
|
||||
@ -60,45 +80,33 @@ Foam::distributionModels::RosinRammler::RosinRammler
|
||||
Foam::distributionModels::RosinRammler::RosinRammler(const RosinRammler& p)
|
||||
:
|
||||
distributionModel(p),
|
||||
minValue_(p.minValue_),
|
||||
maxValue_(p.maxValue_),
|
||||
d_(p.d_),
|
||||
lambda_(p.lambda_),
|
||||
n_(p.n_)
|
||||
{}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
|
||||
|
||||
Foam::distributionModels::RosinRammler::~RosinRammler()
|
||||
{}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
|
||||
|
||||
Foam::scalar Foam::distributionModels::RosinRammler::sample() const
|
||||
{
|
||||
const scalar minValueByDPowN = pow(minValue_/d_, n_);
|
||||
const scalar K = 1 - exp(- pow(maxValue_/d_, n_) + minValueByDPowN);
|
||||
const scalar y = rndGen_.sample01<scalar>();
|
||||
return d_*pow(minValueByDPowN - log(1 - K*y), 1/n_);
|
||||
}
|
||||
|
||||
|
||||
Foam::scalar Foam::distributionModels::RosinRammler::minValue() const
|
||||
{
|
||||
return minValue_;
|
||||
}
|
||||
|
||||
|
||||
Foam::scalar Foam::distributionModels::RosinRammler::maxValue() const
|
||||
{
|
||||
return maxValue_;
|
||||
const scalar u = rndGen_.sample01<scalar>();
|
||||
const scalar qMin = pow(minValue_/lambda_, n_);
|
||||
const scalar qMax = pow(maxValue_/lambda_, n_);
|
||||
const scalar r = scalar(1) - exp(-qMax + qMin);
|
||||
return lambda_*pow(qMin - log(scalar(1) - u*r), scalar(1)/n_);
|
||||
}
|
||||
|
||||
|
||||
Foam::scalar Foam::distributionModels::RosinRammler::meanValue() const
|
||||
{
|
||||
return d_;
|
||||
// (C:Eq. 5)
|
||||
const scalar a = scalar(1)/lambda_ + scalar(1);
|
||||
const scalar qMax = pow(maxValue_/n_, lambda_);
|
||||
const scalar qMin = pow(minValue_/n_, lambda_);
|
||||
const scalar gMax = Math::incGamma_P(a, qMax);
|
||||
const scalar gMin = Math::incGamma_P(a, qMin);
|
||||
|
||||
return n_/(exp(-qMin) - exp(-qMax))*(gMax - gMin);
|
||||
}
|
||||
|
||||
|
||||
|
||||
@ -6,6 +6,7 @@
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
Copyright (C) 2011-2020 OpenFOAM Foundation
|
||||
Copyright (C) 2021 OpenCFD Ltd.
|
||||
-------------------------------------------------------------------------------
|
||||
License
|
||||
This file is part of OpenFOAM.
|
||||
@ -27,13 +28,107 @@ Class
|
||||
Foam::distributionModels::RosinRammler
|
||||
|
||||
Description
|
||||
Rosin-Rammler distributionModel
|
||||
Particle-size distribution model wherein random samples are drawn
|
||||
from the doubly-truncated two-parameter Rosin-Rammler (Weibull)
|
||||
probability density function:
|
||||
|
||||
\f[
|
||||
CDF(x) =
|
||||
(1 - exp(-(x/d)^n + (d_0/d)^n)
|
||||
/(1 - exp(-(d_1/d)^n + (d_0/d)^n)
|
||||
\f]
|
||||
\f[
|
||||
f(x; \lambda, n, A, B) =
|
||||
\frac{
|
||||
\frac{n}{\lambda}
|
||||
\left( \frac{x}{\lambda} \right)^{n-1}
|
||||
\exp\{ -(\frac{x}{\lambda})^n \}
|
||||
}{
|
||||
\exp\{- (\frac{A}{\lambda})^n \}
|
||||
- \exp\{- (\frac{B}{\lambda})^n \}
|
||||
}
|
||||
\f]
|
||||
where
|
||||
|
||||
\vartable
|
||||
f(x; \lambda, n, A, B) | Rosin-Rammler probability density function
|
||||
\lambda | Scale parameter
|
||||
n | Shape parameter
|
||||
x | Sample
|
||||
A | Minimum of the distribution
|
||||
B | Maximum of the distribution
|
||||
\endvartable
|
||||
|
||||
Constraints:
|
||||
- \f$ \infty > B > A > 0\f$
|
||||
- \f$ x \in [B,A] \f$
|
||||
- \f$ \lambda > 0 \f$
|
||||
- \f$ n > 0 \f$
|
||||
|
||||
Random samples are generated by the inverse transform sampling technique
|
||||
by using the quantile function of the doubly-truncated two-parameter
|
||||
Rosin-Rammler (Weibull) probability density function:
|
||||
|
||||
\f[
|
||||
x = \lambda \left( q_{min} - \ln (1 - u r) \right)^{1/n}
|
||||
\f]
|
||||
with
|
||||
|
||||
\f[
|
||||
r = 1 - \exp(-q_{max} + q_{min})
|
||||
\f]
|
||||
|
||||
\f[
|
||||
q_{min} = \left(\frac{A}{\lambda}\right)^n
|
||||
\f]
|
||||
|
||||
\f[
|
||||
q_{max} = \left(\frac{B}{\lambda}\right)^n
|
||||
\f]
|
||||
where \f$ u \f$ is sample drawn from the uniform probability
|
||||
density function on the unit interval \f$ (0, 1) \f$.
|
||||
|
||||
Reference:
|
||||
\verbatim
|
||||
Doubly-truncated two-parameter Weibull distribution and its moments (tag:C):
|
||||
Crénin, F. (2015).
|
||||
Truncated Weibull Distribution Functions and Moments.
|
||||
SSRN 2690255.
|
||||
DOI:10.2139/ssrn.2690255
|
||||
\endverbatim
|
||||
|
||||
Usage
|
||||
Minimal example by using \c constant/\<CloudProperties\>:
|
||||
\verbatim
|
||||
subModels
|
||||
{
|
||||
injectionModels
|
||||
{
|
||||
<name>
|
||||
{
|
||||
...
|
||||
|
||||
sizeDistribution
|
||||
{
|
||||
type RosinRammler;
|
||||
RosinRammlerDistribution
|
||||
{
|
||||
lambda <scaleParameterValue>;
|
||||
n <shapeParameterValue>;
|
||||
minValue <minValue>;
|
||||
maxValue <maxValue>;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
\endverbatim
|
||||
|
||||
where the entries mean:
|
||||
\table
|
||||
Property | Description | Type | Reqd | Deflt
|
||||
type | Type name: RosinRammler | word | yes | -
|
||||
RosinRammlerDistribution | Distribution settings | dict | yes | -
|
||||
lambda | Scale parameter | scalar | yes | -
|
||||
n | Shape parameter | scalar | yes | -
|
||||
minValue | Minimum of the distribution | scalar | yes | -
|
||||
maxValue | Maximum of the distribution | scalar | yes | -
|
||||
\endtable
|
||||
|
||||
SourceFiles
|
||||
RosinRammler.C
|
||||
@ -62,19 +157,11 @@ class RosinRammler
|
||||
{
|
||||
// Private Data
|
||||
|
||||
//- Distribution minimum
|
||||
scalar minValue_;
|
||||
//- Scale parameter
|
||||
scalar lambda_;
|
||||
|
||||
//- Distribution maximum
|
||||
scalar maxValue_;
|
||||
|
||||
// Model coefficients
|
||||
|
||||
//- Scale parameter
|
||||
scalar d_;
|
||||
|
||||
//- Shape parameter
|
||||
scalar n_;
|
||||
//- Shape parameter
|
||||
scalar n_;
|
||||
|
||||
|
||||
public:
|
||||
@ -88,7 +175,7 @@ public:
|
||||
//- Construct from components
|
||||
RosinRammler(const dictionary& dict, Random& rndGen);
|
||||
|
||||
//- Construct copy
|
||||
//- Copy construct
|
||||
RosinRammler(const RosinRammler& p);
|
||||
|
||||
//- Construct and return a clone
|
||||
@ -97,23 +184,20 @@ public:
|
||||
return autoPtr<distributionModel>(new RosinRammler(*this));
|
||||
}
|
||||
|
||||
//- No copy assignment
|
||||
void operator=(const RosinRammler&) = delete;
|
||||
|
||||
|
||||
//- Destructor
|
||||
virtual ~RosinRammler();
|
||||
virtual ~RosinRammler() = default;
|
||||
|
||||
|
||||
// Member Functions
|
||||
|
||||
//- Sample the distributionModel
|
||||
//- Sample the distribution
|
||||
virtual scalar sample() const;
|
||||
|
||||
//- Return the minimum value
|
||||
virtual scalar minValue() const;
|
||||
|
||||
//- Return the maximum value
|
||||
virtual scalar maxValue() const;
|
||||
|
||||
//- Return the mean value
|
||||
//- Return the theoretical mean of the distribution
|
||||
virtual scalar meanValue() const;
|
||||
};
|
||||
|
||||
|
||||
@ -5,7 +5,7 @@
|
||||
\\ / A nd | www.openfoam.com
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
Copyright (C) 2015-2016 OpenCFD Ltd.
|
||||
Copyright (C) 2015-2021 OpenCFD Ltd.
|
||||
-------------------------------------------------------------------------------
|
||||
License
|
||||
This file is part of OpenFOAM.
|
||||
@ -32,11 +32,11 @@ License
|
||||
|
||||
namespace Foam
|
||||
{
|
||||
namespace distributionModels
|
||||
{
|
||||
defineTypeNameAndDebug(binned, 0);
|
||||
addToRunTimeSelectionTable(distributionModel, binned, dictionary);
|
||||
}
|
||||
namespace distributionModels
|
||||
{
|
||||
defineTypeNameAndDebug(binned, 0);
|
||||
addToRunTimeSelectionTable(distributionModel, binned, dictionary);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
@ -58,6 +58,16 @@ void Foam::distributionModels::binned::initialise()
|
||||
|
||||
// Normalise
|
||||
scalar sumProb = xy_.last()[1];
|
||||
|
||||
if (sumProb < VSMALL)
|
||||
{
|
||||
FatalErrorInFunction
|
||||
<< type() << "distribution: "
|
||||
<< "The sum of elements in the second column cannot be zero." << nl
|
||||
<< "sum = " << sumProb
|
||||
<< exit(FatalError);
|
||||
}
|
||||
|
||||
forAll(xy_, bini)
|
||||
{
|
||||
xy_[bini][1] /= sumProb;
|
||||
@ -90,14 +100,10 @@ Foam::distributionModels::binned::binned
|
||||
xy_(distributionModelDict_.lookup("distribution")),
|
||||
meanValue_(0)
|
||||
{
|
||||
if (maxValue() < minValue())
|
||||
{
|
||||
FatalErrorInFunction
|
||||
<< "Maximum value is smaller than the minimum value:"
|
||||
<< " maxValue = " << maxValue()
|
||||
<< ", minValue = " << minValue()
|
||||
<< exit(FatalError);
|
||||
}
|
||||
minValue_ = xy_[0][0];
|
||||
maxValue_ = xy_[xy_.size()-1][0];
|
||||
|
||||
check();
|
||||
|
||||
initialise();
|
||||
}
|
||||
@ -114,16 +120,16 @@ Foam::distributionModels::binned::binned
|
||||
xy_(),
|
||||
meanValue_(0)
|
||||
{
|
||||
scalar minValue = GREAT;
|
||||
scalar maxValue = -GREAT;
|
||||
minValue_ = GREAT;
|
||||
maxValue_ = -GREAT;
|
||||
forAll(sampleData, i)
|
||||
{
|
||||
minValue = min(minValue, sampleData[i]);
|
||||
maxValue = max(maxValue, sampleData[i]);
|
||||
minValue_ = min(minValue_, sampleData[i]);
|
||||
maxValue_ = max(maxValue_, sampleData[i]);
|
||||
}
|
||||
|
||||
const label bin0 = floor(minValue/binWidth);
|
||||
const label bin1 = ceil(maxValue/binWidth);
|
||||
const label bin0 = floor(minValue_/binWidth);
|
||||
const label bin1 = ceil(maxValue_/binWidth);
|
||||
const label nBin = bin1 - bin0;
|
||||
|
||||
if (nBin == 0)
|
||||
@ -177,39 +183,21 @@ Foam::distributionModels::binned::binned(const binned& p)
|
||||
{}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
|
||||
|
||||
Foam::distributionModels::binned::~binned()
|
||||
{}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
|
||||
|
||||
Foam::scalar Foam::distributionModels::binned::sample() const
|
||||
{
|
||||
scalar y = rndGen_.sample01<scalar>();
|
||||
const scalar u = rndGen_.sample01<scalar>();
|
||||
|
||||
for (label i = 0; i < xy_.size() - 1; ++i)
|
||||
{
|
||||
if (xy_[i][1] > y)
|
||||
if (xy_[i][1] > u)
|
||||
{
|
||||
return xy_[i][0];
|
||||
}
|
||||
}
|
||||
|
||||
return xy_.last()[0];
|
||||
}
|
||||
|
||||
|
||||
Foam::scalar Foam::distributionModels::binned::minValue() const
|
||||
{
|
||||
return xy_.first()[0];
|
||||
}
|
||||
|
||||
|
||||
Foam::scalar Foam::distributionModels::binned::maxValue() const
|
||||
{
|
||||
return xy_.last()[0];
|
||||
return maxValue_;
|
||||
}
|
||||
|
||||
|
||||
|
||||
@ -5,7 +5,7 @@
|
||||
\\ / A nd | www.openfoam.com
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
Copyright (C) 2015-2016 OpenCFD Ltd.
|
||||
Copyright (C) 2015-2021 OpenCFD Ltd.
|
||||
-------------------------------------------------------------------------------
|
||||
License
|
||||
This file is part of OpenFOAM.
|
||||
@ -27,8 +27,57 @@ Class
|
||||
Foam::distributionModels::binned
|
||||
|
||||
Description
|
||||
A binned distribution model where the random sample is taken from a
|
||||
discrete set of bins
|
||||
Particle-size distribution model wherein random samples are
|
||||
drawn from a given discrete set of (\c bin, \c probability) pairs,
|
||||
where in terms of its meaning, \c bins correspond to particle sizes
|
||||
and \c probabilities correspond to (relative) probability of occurrences.
|
||||
|
||||
The second column (i.e. \c probability) are normalised by the sum of all
|
||||
its values, resulting in a normalised column in the domain of [0,1].
|
||||
To generate a sample, first a sample drawn from the uniform probability
|
||||
density function on the unit interval (i.e. \c u), and then, the \c bin
|
||||
corresponding to the first \c probability larger than \c u is fetched
|
||||
as the particle size to be further processed.
|
||||
|
||||
Usage
|
||||
Minimal example by using \c constant/\<CloudProperties\>:
|
||||
\verbatim
|
||||
subModels
|
||||
{
|
||||
injectionModels
|
||||
{
|
||||
<name>
|
||||
{
|
||||
...
|
||||
|
||||
sizeDistribution
|
||||
{
|
||||
type binned;
|
||||
binnedDistribution
|
||||
{
|
||||
distribution
|
||||
(
|
||||
(<bin1> <probability1>)
|
||||
(<bin2> <probability2>)
|
||||
...
|
||||
(<binN> <probabilityN>)
|
||||
);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
\endverbatim
|
||||
|
||||
where the entries mean:
|
||||
\table
|
||||
Property | Description | Type | Reqd | Deflt
|
||||
type | Type name: binned | word | yes | -
|
||||
binnedDistribution | Distribution settings | dict | yes | -
|
||||
distribution | \<bin\>-\<probability\> pairs | dict | yes | -
|
||||
\<bin\> | Particle size | scalar | yes | -
|
||||
\<probability\> | Probability of occurrence | scalar | yes | -
|
||||
\endtable
|
||||
|
||||
SourceFiles
|
||||
binned.C
|
||||
@ -70,18 +119,18 @@ class binned
|
||||
:
|
||||
public distributionModel
|
||||
{
|
||||
// Private data
|
||||
typedef VectorSpace<Vector<scalar>, scalar, 2> pair;
|
||||
|
||||
typedef VectorSpace<Vector<scalar>, scalar, 2> pair;
|
||||
// Private Data
|
||||
|
||||
// List of (bin probability)
|
||||
// List of (bin probability) pairs
|
||||
List<pair> xy_;
|
||||
|
||||
//- Distribution mean value
|
||||
//- Mean of the distribution
|
||||
scalar meanValue_;
|
||||
|
||||
|
||||
// Private member functions
|
||||
// Private Member Functions
|
||||
|
||||
//- Initialise the distribution parameters
|
||||
void initialise();
|
||||
@ -101,6 +150,7 @@ public:
|
||||
binned(const dictionary& dict, Random& rndGen);
|
||||
|
||||
//- Construct from components
|
||||
// Allows negative entries
|
||||
binned
|
||||
(
|
||||
const UList<scalar>& sampleData,
|
||||
@ -108,7 +158,7 @@ public:
|
||||
Random& rndGen
|
||||
);
|
||||
|
||||
//- Construct copy
|
||||
//- Copy construct
|
||||
binned(const binned& p);
|
||||
|
||||
//- Construct and return a clone
|
||||
@ -117,23 +167,20 @@ public:
|
||||
return autoPtr<distributionModel>(new binned(*this));
|
||||
}
|
||||
|
||||
//- No copy assignment
|
||||
void operator=(const binned&) = delete;
|
||||
|
||||
|
||||
//- Destructor
|
||||
virtual ~binned();
|
||||
virtual ~binned() = default;
|
||||
|
||||
|
||||
// Member Functions
|
||||
|
||||
//- Sample the distributionModel
|
||||
//- Sample the distribution
|
||||
virtual scalar sample() const;
|
||||
|
||||
//- Return the minimum value
|
||||
virtual scalar minValue() const;
|
||||
|
||||
//- Return the maximum value
|
||||
virtual scalar maxValue() const;
|
||||
|
||||
//- Return the mean value
|
||||
//- Return the arithmetic mean of the distribution data
|
||||
virtual scalar meanValue() const;
|
||||
|
||||
//- Write data to stream
|
||||
|
||||
@ -6,6 +6,7 @@
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
Copyright (C) 2011-2017 OpenFOAM Foundation
|
||||
Copyright (C) 2021 OpenCFD Ltd.
|
||||
-------------------------------------------------------------------------------
|
||||
License
|
||||
This file is part of OpenFOAM.
|
||||
@ -43,19 +44,31 @@ void Foam::distributionModel::check() const
|
||||
if (minValue() < 0)
|
||||
{
|
||||
FatalErrorInFunction
|
||||
<< type() << "distribution: Minimum value must be greater than "
|
||||
<< "zero." << nl << "Supplied minValue = " << minValue()
|
||||
<< type() << "Distribution: "
|
||||
<< "Minimum value must be greater than zero." << nl
|
||||
<< "Supplied minValue = " << minValue()
|
||||
<< abort(FatalError);
|
||||
}
|
||||
|
||||
if (maxValue() < minValue())
|
||||
{
|
||||
FatalErrorInFunction
|
||||
<< type() << "distribution: Maximum value is smaller than the "
|
||||
<< "minimum value:" << nl << " maxValue = " << maxValue()
|
||||
<< type() << "Distribution: "
|
||||
<< "Maximum value cannot be smaller than minimum value" << nl
|
||||
<< " maxValue = " << maxValue()
|
||||
<< ", minValue = " << minValue()
|
||||
<< abort(FatalError);
|
||||
}
|
||||
|
||||
if (maxValue() == minValue())
|
||||
{
|
||||
WarningInFunction
|
||||
<< type() << "Distribution: "
|
||||
<< "Maximum and minimum values are equal to each other" << nl
|
||||
<< " maxValue = " << maxValue()
|
||||
<< ", minValue = " << minValue()
|
||||
<< endl;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
@ -69,7 +82,9 @@ Foam::distributionModel::distributionModel
|
||||
)
|
||||
:
|
||||
distributionModelDict_(dict),
|
||||
rndGen_(rndGen)
|
||||
rndGen_(rndGen),
|
||||
minValue_(distributionModelDict_.getOrDefault<scalar>("minValue", GREAT)),
|
||||
maxValue_(distributionModelDict_.getOrDefault<scalar>("maxValue", -GREAT))
|
||||
{}
|
||||
|
||||
|
||||
@ -79,14 +94,24 @@ Foam::distributionModel::distributionModel
|
||||
)
|
||||
:
|
||||
distributionModelDict_(p.distributionModelDict_),
|
||||
rndGen_(p.rndGen_)
|
||||
rndGen_(p.rndGen_),
|
||||
minValue_(p.minValue_),
|
||||
maxValue_(p.maxValue_)
|
||||
{}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
|
||||
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
|
||||
|
||||
Foam::distributionModel::~distributionModel()
|
||||
{}
|
||||
Foam::scalar Foam::distributionModel::minValue() const
|
||||
{
|
||||
return minValue_;
|
||||
}
|
||||
|
||||
|
||||
Foam::scalar Foam::distributionModel::maxValue() const
|
||||
{
|
||||
return maxValue_;
|
||||
}
|
||||
|
||||
|
||||
// ************************************************************************* //
|
||||
|
||||
@ -6,6 +6,7 @@
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
Copyright (C) 2011-2017 OpenFOAM Foundation
|
||||
Copyright (C) 2021 OpenCFD Ltd.
|
||||
-------------------------------------------------------------------------------
|
||||
License
|
||||
This file is part of OpenFOAM.
|
||||
@ -23,15 +24,21 @@ License
|
||||
You should have received a copy of the GNU General Public License
|
||||
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
|
||||
|
||||
Namespace
|
||||
Foam::distributionModels
|
||||
|
||||
Description
|
||||
A namespace for various probability distribution model implementations.
|
||||
|
||||
Class
|
||||
Foam::distributionModel
|
||||
|
||||
Description
|
||||
A library of runtime-selectable distribution models.
|
||||
A library of runtime-selectable doubly-truncated probability distribution
|
||||
models. Returns random samples based on given distribution parameters.
|
||||
|
||||
Returns a sampled value given the expectation (nu) and variance (sigma^2)
|
||||
|
||||
Current distribution models include:
|
||||
Available distribution models include:
|
||||
- binned
|
||||
- exponential
|
||||
- fixedValue
|
||||
- general
|
||||
@ -41,11 +48,6 @@ Description
|
||||
- Mass-based Rosin-Rammler
|
||||
- uniform
|
||||
|
||||
The distributionModel is tabulated in equidistant nPoints, in an interval.
|
||||
These values are integrated to obtain the cumulated distribution model,
|
||||
which is then used to change the distribution from unifrom to
|
||||
the actual distributionModel.
|
||||
|
||||
SourceFiles
|
||||
distributionModel.C
|
||||
distributionModelNew.C
|
||||
@ -70,10 +72,9 @@ namespace Foam
|
||||
|
||||
class distributionModel
|
||||
{
|
||||
|
||||
protected:
|
||||
|
||||
// Protected data
|
||||
// Protected Data
|
||||
|
||||
//- Coefficients dictionary
|
||||
const dictionary distributionModelDict_;
|
||||
@ -81,6 +82,12 @@ protected:
|
||||
//- Reference to the random number generator
|
||||
Random& rndGen_;
|
||||
|
||||
//- Minimum of the distribution
|
||||
scalar minValue_;
|
||||
|
||||
//- Maximum of the distribution
|
||||
scalar maxValue_;
|
||||
|
||||
|
||||
// Protected Member Functions
|
||||
|
||||
@ -118,7 +125,7 @@ public:
|
||||
Random& rndGen
|
||||
);
|
||||
|
||||
//- Construct copy
|
||||
//- Copy construct
|
||||
distributionModel(const distributionModel& p);
|
||||
|
||||
//- Construct and return a clone
|
||||
@ -134,21 +141,22 @@ public:
|
||||
|
||||
|
||||
//- Destructor
|
||||
virtual ~distributionModel();
|
||||
virtual ~distributionModel() = default;
|
||||
|
||||
|
||||
// Member Functions
|
||||
|
||||
//- Sample the distributionModel
|
||||
//- Sample the distribution
|
||||
virtual scalar sample() const = 0;
|
||||
|
||||
//- Return the minimum value
|
||||
virtual scalar minValue() const = 0;
|
||||
//- Return the minimum of the distribution
|
||||
virtual scalar minValue() const;
|
||||
|
||||
//- Return the maximum value
|
||||
virtual scalar maxValue() const = 0;
|
||||
//- Return the maximum of the distribution
|
||||
virtual scalar maxValue() const;
|
||||
|
||||
//- Return the maximum value
|
||||
//- Return the theoretical mean of the distribution, or
|
||||
//- the arithmetic mean of the distribution data
|
||||
virtual scalar meanValue() const = 0;
|
||||
};
|
||||
|
||||
|
||||
@ -6,6 +6,7 @@
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
Copyright (C) 2011-2016 OpenFOAM Foundation
|
||||
Copyright (C) 2021 OpenCFD Ltd.
|
||||
-------------------------------------------------------------------------------
|
||||
License
|
||||
This file is part of OpenFOAM.
|
||||
@ -48,10 +49,16 @@ Foam::distributionModels::exponential::exponential
|
||||
)
|
||||
:
|
||||
distributionModel(typeName, dict, rndGen),
|
||||
minValue_(distributionModelDict_.get<scalar>("minValue")),
|
||||
maxValue_(distributionModelDict_.get<scalar>("maxValue")),
|
||||
lambda_(distributionModelDict_.get<scalar>("lambda"))
|
||||
{
|
||||
if (lambda_ < VSMALL)
|
||||
{
|
||||
FatalErrorInFunction
|
||||
<< "Rate parameter cannot be equal to or less than zero:" << nl
|
||||
<< " lambda = " << lambda_
|
||||
<< exit(FatalError);
|
||||
}
|
||||
|
||||
check();
|
||||
}
|
||||
|
||||
@ -59,43 +66,24 @@ Foam::distributionModels::exponential::exponential
|
||||
Foam::distributionModels::exponential::exponential(const exponential& p)
|
||||
:
|
||||
distributionModel(p),
|
||||
minValue_(p.minValue_),
|
||||
maxValue_(p.maxValue_),
|
||||
lambda_(p.lambda_)
|
||||
{}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
|
||||
|
||||
Foam::distributionModels::exponential::~exponential()
|
||||
{}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
|
||||
|
||||
Foam::scalar Foam::distributionModels::exponential::sample() const
|
||||
{
|
||||
scalar y = rndGen_.sample01<scalar>();
|
||||
scalar K = exp(-lambda_*maxValue_) - exp(-lambda_*minValue_);
|
||||
return -(1.0/lambda_)*log(exp(-lambda_*minValue_) + y*K);
|
||||
}
|
||||
|
||||
|
||||
Foam::scalar Foam::distributionModels::exponential::minValue() const
|
||||
{
|
||||
return minValue_;
|
||||
}
|
||||
|
||||
|
||||
Foam::scalar Foam::distributionModels::exponential::maxValue() const
|
||||
{
|
||||
return maxValue_;
|
||||
const scalar u = rndGen_.sample01<scalar>();
|
||||
const scalar qMin = exp(-lambda_*minValue_);
|
||||
const scalar qMax = exp(-lambda_*maxValue_);
|
||||
return -(scalar(1)/lambda_)*log(qMin + u*(qMax - qMin));
|
||||
}
|
||||
|
||||
|
||||
Foam::scalar Foam::distributionModels::exponential::meanValue() const
|
||||
{
|
||||
return 1.0/lambda_;
|
||||
return scalar(1)/lambda_;
|
||||
}
|
||||
|
||||
|
||||
|
||||
@ -6,6 +6,7 @@
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
Copyright (C) 2011-2016 OpenFOAM Foundation
|
||||
Copyright (C) 2021 OpenCFD Ltd.
|
||||
-------------------------------------------------------------------------------
|
||||
License
|
||||
This file is part of OpenFOAM.
|
||||
@ -27,7 +28,87 @@ Class
|
||||
Foam::distributionModels::exponential
|
||||
|
||||
Description
|
||||
exponential distribution model
|
||||
Particle-size distribution model wherein random samples are drawn
|
||||
from the doubly-truncated exponential probability density function:
|
||||
|
||||
\f[
|
||||
f(x; \lambda, A, B) =
|
||||
\lambda \frac{\exp(-\lambda (x - A))}{1 - \exp(-\lambda(B-A))}
|
||||
\f]
|
||||
where
|
||||
|
||||
\vartable
|
||||
f(x; \lambda, A, B) | Exponential probability density function
|
||||
\lambda | Rate parameter
|
||||
x | Sample
|
||||
A | Minimum of the distribution
|
||||
B | Maximum of the distribution
|
||||
\endvartable
|
||||
|
||||
Constraints:
|
||||
- \f$ \infty > B > A > 0 \f$
|
||||
- \f$ x \in [B,A] \f$
|
||||
- \f$ \lambda > 0 \f$
|
||||
|
||||
Random samples are generated by the inverse transform sampling technique
|
||||
by using the quantile function of the doubly-truncated exponential
|
||||
probability density function:
|
||||
|
||||
\f[
|
||||
x = - \frac{1}{\lambda} \ln (r)
|
||||
\f]
|
||||
|
||||
with
|
||||
|
||||
\f[
|
||||
r = q_{min} + u (q_{max} - q_{min})
|
||||
\f]
|
||||
|
||||
\f[
|
||||
q_{min} = \exp(-\lambda A)
|
||||
\f]
|
||||
|
||||
\f[
|
||||
q_{max} = \exp(-\lambda B)
|
||||
\f]
|
||||
where \f$ u \f$ is a sample drawn from the uniform probability
|
||||
density function on the unit interval \f$ (0, 1) \f$.
|
||||
|
||||
Usage
|
||||
Minimal example by using \c constant/\<CloudProperties\>:
|
||||
\verbatim
|
||||
subModels
|
||||
{
|
||||
injectionModels
|
||||
{
|
||||
<name>
|
||||
{
|
||||
...
|
||||
|
||||
sizeDistribution
|
||||
{
|
||||
type exponential;
|
||||
exponentialDistribution
|
||||
{
|
||||
lambda <lambdaValue>;
|
||||
minValue <minValue>;
|
||||
maxValue <maxValue>;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
\endverbatim
|
||||
|
||||
where the entries mean:
|
||||
\table
|
||||
Property | Description | Type | Reqd | Deflt
|
||||
type | Type name: exponential | word | yes | -
|
||||
exponentialDistribution | Distribution settings | dict | yes | -
|
||||
lambda | Rate parameter | scalar | yes | -
|
||||
minValue | Minimum of the distribution | scalar | yes | -
|
||||
maxValue | Maximum of the distribution | scalar | yes | -
|
||||
\endtable
|
||||
|
||||
SourceFiles
|
||||
exponential.C
|
||||
@ -54,18 +135,10 @@ class exponential
|
||||
:
|
||||
public distributionModel
|
||||
{
|
||||
// Private data
|
||||
// Private Data
|
||||
|
||||
//- Distribution minimum
|
||||
scalar minValue_;
|
||||
|
||||
//- Distribution maximum
|
||||
scalar maxValue_;
|
||||
|
||||
|
||||
// Model coefficients
|
||||
|
||||
scalar lambda_;
|
||||
//- Rate parameter
|
||||
scalar lambda_;
|
||||
|
||||
|
||||
public:
|
||||
@ -79,7 +152,7 @@ public:
|
||||
//- Construct from components
|
||||
exponential(const dictionary& dict, Random& rndGen);
|
||||
|
||||
//- Construct copy
|
||||
//- Copy construct
|
||||
exponential(const exponential& p);
|
||||
|
||||
//- Construct and return a clone
|
||||
@ -88,23 +161,20 @@ public:
|
||||
return autoPtr<distributionModel>(new exponential(*this));
|
||||
}
|
||||
|
||||
//- No copy assignment
|
||||
void operator=(const exponential&) = delete;
|
||||
|
||||
|
||||
//- Destructor
|
||||
virtual ~exponential();
|
||||
virtual ~exponential() = default;
|
||||
|
||||
|
||||
// Member Functions
|
||||
|
||||
//- Sample the distributionModel
|
||||
//- Sample the distribution
|
||||
virtual scalar sample() const;
|
||||
|
||||
//- Return the minimum value
|
||||
virtual scalar minValue() const;
|
||||
|
||||
//- Return the maximum value
|
||||
virtual scalar maxValue() const;
|
||||
|
||||
//- Return the mean value
|
||||
//- Return the theoretical mean of the distribution
|
||||
virtual scalar meanValue() const;
|
||||
};
|
||||
|
||||
|
||||
@ -6,6 +6,7 @@
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
Copyright (C) 2011-2016 OpenFOAM Foundation
|
||||
Copyright (C) 2021 OpenCFD Ltd.
|
||||
-------------------------------------------------------------------------------
|
||||
License
|
||||
This file is part of OpenFOAM.
|
||||
@ -49,7 +50,15 @@ Foam::distributionModels::fixedValue::fixedValue
|
||||
:
|
||||
distributionModel(typeName, dict, rndGen),
|
||||
value_(distributionModelDict_.get<scalar>("value"))
|
||||
{}
|
||||
{
|
||||
if (value_ < VSMALL)
|
||||
{
|
||||
FatalErrorInFunction
|
||||
<< "Fixed value cannot be equal to or less than zero:" << nl
|
||||
<< " value = " << value_
|
||||
<< exit(FatalError);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
Foam::distributionModels::fixedValue::fixedValue(const fixedValue& p)
|
||||
@ -59,12 +68,6 @@ Foam::distributionModels::fixedValue::fixedValue(const fixedValue& p)
|
||||
{}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
|
||||
|
||||
Foam::distributionModels::fixedValue::~fixedValue()
|
||||
{}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
|
||||
|
||||
Foam::scalar Foam::distributionModels::fixedValue::fixedValue::sample() const
|
||||
|
||||
@ -6,6 +6,7 @@
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
Copyright (C) 2011-2016 OpenFOAM Foundation
|
||||
Copyright (C) 2021 OpenCFD Ltd.
|
||||
-------------------------------------------------------------------------------
|
||||
License
|
||||
This file is part of OpenFOAM.
|
||||
@ -27,7 +28,39 @@ Class
|
||||
Foam::distributionModels::fixedValue
|
||||
|
||||
Description
|
||||
Returns a fixed value
|
||||
Particle-size distribution model wherein samples are given fixed values.
|
||||
|
||||
Usage
|
||||
Minimal example by using \c constant/\<CloudProperties\>:
|
||||
\verbatim
|
||||
subModels
|
||||
{
|
||||
injectionModels
|
||||
{
|
||||
<name>
|
||||
{
|
||||
...
|
||||
|
||||
sizeDistribution
|
||||
{
|
||||
type fixedValue;
|
||||
fixedValueDistribution
|
||||
{
|
||||
value <value>;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
\endverbatim
|
||||
|
||||
where the entries mean:
|
||||
\table
|
||||
Property | Description | Type | Reqd | Deflt
|
||||
type | Type name: fixedValue | word | yes | -
|
||||
fixedValueDistribution | Distribution settings | dict | yes | -
|
||||
value | Fixed value for size | scalar | yes | -
|
||||
\endtable
|
||||
|
||||
SourceFiles
|
||||
fixedValue.C
|
||||
@ -53,9 +86,9 @@ class fixedValue
|
||||
:
|
||||
public distributionModel
|
||||
{
|
||||
// Private data
|
||||
// Private Data
|
||||
|
||||
//- Fixed value
|
||||
//- Fixed value for size
|
||||
scalar value_;
|
||||
|
||||
|
||||
@ -70,7 +103,7 @@ public:
|
||||
//- Construct from components
|
||||
fixedValue(const dictionary& dict, Random& rndGen);
|
||||
|
||||
//- Construct copy
|
||||
//- Copy construct
|
||||
fixedValue(const fixedValue& p);
|
||||
|
||||
//- Construct and return a clone
|
||||
@ -79,23 +112,26 @@ public:
|
||||
return autoPtr<distributionModel>(new fixedValue(*this));
|
||||
}
|
||||
|
||||
//- No copy assignment
|
||||
void operator=(const fixedValue&) = delete;
|
||||
|
||||
|
||||
//- Destructor
|
||||
virtual ~fixedValue();
|
||||
virtual ~fixedValue() = default;
|
||||
|
||||
|
||||
// Member Functions
|
||||
|
||||
//- Sample the distributionModel
|
||||
//- Sample the distribution
|
||||
virtual scalar sample() const;
|
||||
|
||||
//- Return the minimum value
|
||||
//- Return the minimum of the distribution
|
||||
virtual scalar minValue() const;
|
||||
|
||||
//- Return the maximum value
|
||||
//- Return the maximum of the distribution
|
||||
virtual scalar maxValue() const;
|
||||
|
||||
//- Return the mean value
|
||||
//- Return the theoretical mean of the distribution
|
||||
virtual scalar meanValue() const;
|
||||
};
|
||||
|
||||
|
||||
@ -5,8 +5,8 @@
|
||||
\\ / A nd | www.openfoam.com
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
Copyright (C) 2011-2016 OpenFOAM Foundation
|
||||
Copyright (C) 2016 OpenCFD Ltd.
|
||||
Copyright (C) 2011-2019 OpenFOAM Foundation
|
||||
Copyright (C) 2016-2021 OpenCFD Ltd.
|
||||
-------------------------------------------------------------------------------
|
||||
License
|
||||
This file is part of OpenFOAM.
|
||||
@ -47,13 +47,14 @@ void Foam::distributionModels::general::initialise()
|
||||
{
|
||||
static scalar eps = ROOTVSMALL;
|
||||
|
||||
const label nEntries = xy_.size();
|
||||
integral_.setSize(nEntries_);
|
||||
|
||||
integral_.setSize(nEntries);
|
||||
// Fill out the integral table (x, P(x<=0)) and calculate mean
|
||||
// For density function: P(x<=0) = int f(x) and mean = int x*f(x)
|
||||
// For cumulative function: mean = int 1-P(x<=0) = maxValue_ - int P(x<=0)
|
||||
integral_[0] = 0;
|
||||
|
||||
// Normalize the cumulative distribution
|
||||
integral_[0] = 0.0;
|
||||
for (label i = 1; i < nEntries; i++)
|
||||
for (label i = 1; i < nEntries_; ++i)
|
||||
{
|
||||
scalar k = (xy_[i][1] - xy_[i-1][1])/(xy_[i][0] - xy_[i-1][0] + eps);
|
||||
scalar d = xy_[i-1][1] - k*xy_[i-1][0];
|
||||
@ -61,18 +62,32 @@ void Foam::distributionModels::general::initialise()
|
||||
scalar y0 = xy_[i-1][0]*(0.5*k*xy_[i-1][0] + d);
|
||||
scalar area = y1 - y0;
|
||||
|
||||
integral_[i] = area + integral_[i-1];
|
||||
if (cumulative_)
|
||||
{
|
||||
integral_[i] = xy_[i][1];
|
||||
meanValue_ += area;
|
||||
}
|
||||
else
|
||||
{
|
||||
integral_[i] = area + integral_[i-1];
|
||||
|
||||
y1 = sqr(xy_[i][0])*(1.0/3.0*k*xy_[i][0] + 0.5*d);
|
||||
y0 = sqr(xy_[i-1][0])*(1.0/3.0*k*xy_[i-1][0] + 0.5*d);
|
||||
meanValue_ += y1 - y0;
|
||||
}
|
||||
}
|
||||
|
||||
// normalize the distribution
|
||||
scalar sumArea = integral_.last();
|
||||
|
||||
meanValue_ = sumArea/(maxValue() - minValue() + eps);
|
||||
|
||||
for (label i=0; i < nEntries; i++)
|
||||
for (label i = 0; i < nEntries_; ++i)
|
||||
{
|
||||
xy_[i][1] /= sumArea + eps;
|
||||
integral_[i] /= sumArea + eps;
|
||||
}
|
||||
|
||||
meanValue_ /= sumArea;
|
||||
meanValue_ = cumulative_ ? (maxValue_ - meanValue_ + eps) : meanValue_;
|
||||
}
|
||||
|
||||
|
||||
@ -86,11 +101,64 @@ Foam::distributionModels::general::general
|
||||
:
|
||||
distributionModel(typeName, dict, rndGen),
|
||||
xy_(distributionModelDict_.lookup("distribution")),
|
||||
meanValue_(0.0),
|
||||
integral_()
|
||||
nEntries_(xy_.size()),
|
||||
meanValue_(0),
|
||||
integral_(nEntries_),
|
||||
cumulative_(distributionModelDict_.getOrDefault("cumulative", false))
|
||||
{
|
||||
minValue_ = xy_[0][0];
|
||||
maxValue_ = xy_[nEntries_-1][0];
|
||||
|
||||
check();
|
||||
|
||||
// Additional sanity checks
|
||||
if (cumulative_ && xy_[0][1] != 0)
|
||||
{
|
||||
FatalErrorInFunction
|
||||
<< type() << "distribution: "
|
||||
<< "Elements in the second column for cumulative "
|
||||
<< "distribution functions must start from zero." << nl
|
||||
<< "First element = " << xy_[0][1]
|
||||
<< exit(FatalError);
|
||||
}
|
||||
|
||||
for (label i = 0; i < nEntries_; ++i)
|
||||
{
|
||||
if (i > 0 && xy_[i][0] <= xy_[i-1][0])
|
||||
{
|
||||
FatalErrorInFunction
|
||||
<< type() << "distribution: "
|
||||
<< "Elements in the first column must "
|
||||
<< "be specified in an ascending order." << nl
|
||||
<< "Please see the row i = " << i << nl
|
||||
<< "xy[i-1] = " << xy_[i-1] << nl
|
||||
<< "xy[i] = " << xy_[i]
|
||||
<< exit(FatalError);
|
||||
}
|
||||
|
||||
if (xy_[i][0] < 0 || xy_[i][1] < 0)
|
||||
{
|
||||
FatalErrorInFunction
|
||||
<< type() << "distribution: "
|
||||
<< "Input pairs cannot contain any negative element." << nl
|
||||
<< "Please see the row i = " << i << nl
|
||||
<< "xy[i] = " << xy_[i]
|
||||
<< exit(FatalError);
|
||||
}
|
||||
|
||||
if (cumulative_ && i > 0 && xy_[i][1] < xy_[i-1][1])
|
||||
{
|
||||
FatalErrorInFunction
|
||||
<< type() << "distribution: "
|
||||
<< "Elements in the second column for cumulative "
|
||||
<< "distribution functions must be non-decreasing." << nl
|
||||
<< "Please see the row i = " << i << nl
|
||||
<< "xy[i-1] = " << xy_[i-1] << nl
|
||||
<< "xy[i] = " << xy_[i]
|
||||
<< exit(FatalError);
|
||||
}
|
||||
}
|
||||
|
||||
initialise();
|
||||
}
|
||||
|
||||
@ -104,22 +172,24 @@ Foam::distributionModels::general::general
|
||||
:
|
||||
distributionModel(typeName, dictionary::null, rndGen),
|
||||
xy_(),
|
||||
meanValue_(0.0),
|
||||
integral_()
|
||||
nEntries_(0),
|
||||
meanValue_(0),
|
||||
integral_(),
|
||||
cumulative_(false)
|
||||
{
|
||||
scalar minValue = GREAT;
|
||||
scalar maxValue = -GREAT;
|
||||
minValue_ = GREAT;
|
||||
maxValue_ = -GREAT;
|
||||
forAll(sampleData, i)
|
||||
{
|
||||
minValue = min(minValue, sampleData[i]);
|
||||
maxValue = max(maxValue, sampleData[i]);
|
||||
minValue_ = min(minValue_, sampleData[i]);
|
||||
maxValue_ = max(maxValue_, sampleData[i]);
|
||||
}
|
||||
|
||||
label bin0 = floor(minValue/binWidth);
|
||||
label bin1 = ceil(maxValue/binWidth);
|
||||
label nEntries = bin1 - bin0;
|
||||
label bin0 = floor(minValue_/binWidth);
|
||||
label bin1 = ceil(maxValue_/binWidth);
|
||||
nEntries_ = bin1 - bin0;
|
||||
|
||||
if (nEntries == 0)
|
||||
if (nEntries_ == 0)
|
||||
{
|
||||
WarningInFunction
|
||||
<< "Data cannot be binned - zero bins generated" << nl
|
||||
@ -130,10 +200,10 @@ Foam::distributionModels::general::general
|
||||
return;
|
||||
}
|
||||
|
||||
xy_.setSize(nEntries);
|
||||
xy_.setSize(nEntries_);
|
||||
|
||||
// Populate bin boundaries and initialise occurrences
|
||||
for (label bini = 0; bini < nEntries; ++bini)
|
||||
for (label bini = 0; bini < nEntries_; ++bini)
|
||||
{
|
||||
xy_[bini][0] = (bin0 + bini)*binWidth;
|
||||
xy_[bini][1] = 0;
|
||||
@ -154,14 +224,10 @@ Foam::distributionModels::general::general(const general& p)
|
||||
:
|
||||
distributionModel(p),
|
||||
xy_(p.xy_),
|
||||
nEntries_(p.nEntries_),
|
||||
meanValue_(p.meanValue_),
|
||||
integral_(p.integral_)
|
||||
{}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
|
||||
|
||||
Foam::distributionModels::general::~general()
|
||||
integral_(p.integral_),
|
||||
cumulative_(p.cumulative_)
|
||||
{}
|
||||
|
||||
|
||||
@ -169,57 +235,46 @@ Foam::distributionModels::general::~general()
|
||||
|
||||
Foam::scalar Foam::distributionModels::general::sample() const
|
||||
{
|
||||
scalar y = rndGen_.sample01<scalar>();
|
||||
const scalar u = rndGen_.sample01<scalar>();
|
||||
|
||||
// Find the interval where y is in the table
|
||||
// Find the interval where u is in the table
|
||||
label n = 1;
|
||||
while (integral_[n] <= y)
|
||||
while (integral_[n] <= u)
|
||||
{
|
||||
n++;
|
||||
}
|
||||
|
||||
scalar k = (xy_[n][1] - xy_[n-1][1])/(xy_[n][0] - xy_[n-1][0]);
|
||||
scalar d = xy_[n-1][1] - k*xy_[n-1][0];
|
||||
const scalar k = (xy_[n][1] - xy_[n-1][1])/(xy_[n][0] - xy_[n-1][0]);
|
||||
const scalar d = xy_[n-1][1] - k*xy_[n-1][0];
|
||||
|
||||
scalar alpha = y + xy_[n-1][0]*(0.5*k*xy_[n-1][0] + d) - integral_[n-1];
|
||||
scalar x = 0.0;
|
||||
if (cumulative_)
|
||||
{
|
||||
return (u - d)/k;
|
||||
}
|
||||
|
||||
const scalar alpha =
|
||||
u + xy_[n-1][0]*(0.5*k*xy_[n-1][0] + d) - integral_[n-1];
|
||||
|
||||
// If k is small it is a linear equation, otherwise it is of second order
|
||||
if (mag(k) > SMALL)
|
||||
{
|
||||
scalar p = 2.0*d/k;
|
||||
scalar q = -2.0*alpha/k;
|
||||
scalar sqrtEr = sqrt(0.25*p*p - q);
|
||||
const scalar p = 2.0*d/k;
|
||||
const scalar q = -2.0*alpha/k;
|
||||
const scalar sqrtEr = sqrt(0.25*p*p - q);
|
||||
|
||||
scalar x1 = -0.5*p + sqrtEr;
|
||||
scalar x2 = -0.5*p - sqrtEr;
|
||||
const scalar x1 = -0.5*p + sqrtEr;
|
||||
const scalar x2 = -0.5*p - sqrtEr;
|
||||
if ((x1 >= xy_[n-1][0]) && (x1 <= xy_[n][0]))
|
||||
{
|
||||
x = x1;
|
||||
return x1;
|
||||
}
|
||||
else
|
||||
{
|
||||
x = x2;
|
||||
return x2;
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
x = alpha/d;
|
||||
}
|
||||
|
||||
return x;
|
||||
}
|
||||
|
||||
|
||||
Foam::scalar Foam::distributionModels::general::minValue() const
|
||||
{
|
||||
return xy_.first()[0];
|
||||
}
|
||||
|
||||
|
||||
Foam::scalar Foam::distributionModels::general::maxValue() const
|
||||
{
|
||||
return xy_.last()[0];
|
||||
return alpha/d;
|
||||
}
|
||||
|
||||
|
||||
@ -278,11 +333,11 @@ void Foam::distributionModels::general::readDict(const dictionary& dict)
|
||||
Foam::tmp<Foam::Field<Foam::scalar>>
|
||||
Foam::distributionModels::general::x() const
|
||||
{
|
||||
tmp<Field<scalar>> tx(new Field<scalar>(xy_.size()));
|
||||
scalarField& xi = tx.ref();
|
||||
auto tx = tmp<scalarField>::New(xy_.size());
|
||||
auto& x = tx.ref();
|
||||
forAll(xy_, i)
|
||||
{
|
||||
xi[i] = xy_[i][0];
|
||||
x[i] = xy_[i][0];
|
||||
}
|
||||
|
||||
return tx;
|
||||
@ -292,11 +347,11 @@ Foam::distributionModels::general::x() const
|
||||
Foam::tmp<Foam::Field<Foam::scalar>>
|
||||
Foam::distributionModels::general::y() const
|
||||
{
|
||||
tmp<Field<scalar>> ty(new Field<scalar>(xy_.size()));
|
||||
scalarField& yi = ty.ref();
|
||||
auto ty = tmp<scalarField>::New(xy_.size());
|
||||
auto& y = ty.ref();
|
||||
forAll(xy_, i)
|
||||
{
|
||||
yi[i] = xy_[i][1];
|
||||
y[i] = xy_[i][1];
|
||||
}
|
||||
|
||||
return ty;
|
||||
|
||||
@ -5,8 +5,8 @@
|
||||
\\ / A nd | www.openfoam.com
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
Copyright (C) 2011-2015 OpenFOAM Foundation
|
||||
Copyright (C) 2016 OpenCFD Ltd.
|
||||
Copyright (C) 2011-2019 OpenFOAM Foundation
|
||||
Copyright (C) 2016-2021 OpenCFD Ltd.
|
||||
-------------------------------------------------------------------------------
|
||||
License
|
||||
This file is part of OpenFOAM.
|
||||
@ -28,7 +28,66 @@ Class
|
||||
Foam::distributionModels::general
|
||||
|
||||
Description
|
||||
general distribution model
|
||||
Particle-size distribution model wherein random samples are
|
||||
drawn from a given arbitrary probability density function
|
||||
or cumulative distribution function. Input distributions are
|
||||
specified as pairs of (size, probability) (i.e. (point, value) ).
|
||||
|
||||
Usage
|
||||
Minimal example by using \c constant/\<CloudProperties\>:
|
||||
\verbatim
|
||||
subModels
|
||||
{
|
||||
injectionModels
|
||||
{
|
||||
<name>
|
||||
{
|
||||
...
|
||||
|
||||
sizeDistribution
|
||||
{
|
||||
type general;
|
||||
generalDistribution
|
||||
{
|
||||
cumulative false;
|
||||
|
||||
distribution
|
||||
(
|
||||
(<size1> <probability1>)
|
||||
(<size2> <probability2>)
|
||||
...
|
||||
(<sizeN> <probabilityN>)
|
||||
);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
\endverbatim
|
||||
|
||||
where the entries mean:
|
||||
\table
|
||||
Property | Description | Type | Reqd | Deflt
|
||||
type | Type name: general | word | yes | -
|
||||
generalDistribution | Distribution settings | dict | yes | -
|
||||
distribution | \<size\>-\<probability\> pairs | dict | yes | -
|
||||
\<size\> | Particle size | scalar | yes | -
|
||||
\<probability\> | Volume fraction/probability | scalar | yes | -
|
||||
cumulative | Flag to determine if input distribution is specified <!--
|
||||
--> as cumulative or as density | bool | no | false
|
||||
\endtable
|
||||
|
||||
Note
|
||||
- An example for a pair within \c distribution subdictionary
|
||||
can be given as follows: "(1e-07 1.3)" means 1.3\% of particles
|
||||
are modelled with a particle diameter of "1e-7" [m], and so on.
|
||||
- Variation between input pairs is assumed to be linear.
|
||||
- Elements in the second column (i.e. probability) are normalised.
|
||||
- Elements in the second column for cumulative distribution
|
||||
functions must start from zero and must be non-decreasing (i.e. monotonic).
|
||||
- Elements in the first column (i.e. size) must be specified
|
||||
in an ascending order.
|
||||
- Input pairs cannot contain any negative element.
|
||||
|
||||
SourceFiles
|
||||
general.C
|
||||
@ -65,26 +124,34 @@ namespace distributionModels
|
||||
{
|
||||
|
||||
/*---------------------------------------------------------------------------*\
|
||||
Class general Declaration
|
||||
Class general Declaration
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
class general
|
||||
:
|
||||
public distributionModel
|
||||
{
|
||||
// Private data
|
||||
typedef VectorSpace<Vector<scalar>, scalar, 2> pair;
|
||||
|
||||
typedef VectorSpace<Vector<scalar>, scalar, 2> pair;
|
||||
// Private Data
|
||||
|
||||
// List of (bin probability)
|
||||
//- List of (x, y=f(x)) pairs
|
||||
List<pair> xy_;
|
||||
|
||||
//- Number of entries in the xy_ list
|
||||
label nEntries_;
|
||||
|
||||
//- Mean of the distribution
|
||||
scalar meanValue_;
|
||||
|
||||
//- Values of cumulative distribution function
|
||||
List<scalar> integral_;
|
||||
|
||||
//- Flag to determine if input is specified as cumulative or as density
|
||||
bool cumulative_;
|
||||
|
||||
// Private member functions
|
||||
|
||||
// Private Member Functions
|
||||
|
||||
//- Initialise the distribution parameters
|
||||
void initialise();
|
||||
@ -102,6 +169,7 @@ public:
|
||||
general(const dictionary& dict, Random& rndGen);
|
||||
|
||||
//- Construct from components
|
||||
// Allows negative entries
|
||||
general
|
||||
(
|
||||
const UList<scalar>& sampleData,
|
||||
@ -109,7 +177,7 @@ public:
|
||||
Random& rndGen
|
||||
);
|
||||
|
||||
//- Construct copy
|
||||
//- Copy construct
|
||||
general(const general& p);
|
||||
|
||||
//- Construct and return a clone
|
||||
@ -118,29 +186,26 @@ public:
|
||||
return autoPtr<distributionModel>(new general(*this));
|
||||
}
|
||||
|
||||
//- No copy assignment
|
||||
void operator=(const general&) = delete;
|
||||
|
||||
|
||||
//- Destructor
|
||||
virtual ~general();
|
||||
virtual ~general() = default;
|
||||
|
||||
|
||||
// Member Functions
|
||||
|
||||
//- Bin boundaries
|
||||
virtual tmp<Field<scalar>> x() const;
|
||||
virtual tmp<scalarField> x() const;
|
||||
|
||||
//- Probabilities
|
||||
virtual tmp<Field<scalar>> y() const;
|
||||
virtual tmp<scalarField> y() const;
|
||||
|
||||
//- Sample the distributionModel
|
||||
//- Sample the distribution
|
||||
virtual scalar sample() const;
|
||||
|
||||
//- Return the minimum value
|
||||
virtual scalar minValue() const;
|
||||
|
||||
//- Return the maximum value
|
||||
virtual scalar maxValue() const;
|
||||
|
||||
//- Return the mean value
|
||||
//- Return the arithmetic mean of the distribution data
|
||||
virtual scalar meanValue() const;
|
||||
|
||||
//- Write data to stream
|
||||
|
||||
@ -50,11 +50,18 @@ Foam::distributionModels::massRosinRammler::massRosinRammler
|
||||
)
|
||||
:
|
||||
distributionModel(typeName, dict, rndGen),
|
||||
minValue_(distributionModelDict_.get<scalar>("minValue")),
|
||||
maxValue_(distributionModelDict_.get<scalar>("maxValue")),
|
||||
d_(distributionModelDict_.get<scalar>("d")),
|
||||
lambda_(distributionModelDict_.getCompat<scalar>("lambda", {{"d", 2112}})),
|
||||
n_(distributionModelDict_.get<scalar>("n"))
|
||||
{
|
||||
if (lambda_ < VSMALL || n_ < VSMALL)
|
||||
{
|
||||
FatalErrorInFunction
|
||||
<< "Scale/Shape parameter cannot be equal to or less than zero:"
|
||||
<< " lambda = " << lambda_
|
||||
<< " n = " << n_
|
||||
<< exit(FatalError);
|
||||
}
|
||||
|
||||
check();
|
||||
}
|
||||
|
||||
@ -65,53 +72,37 @@ Foam::distributionModels::massRosinRammler::massRosinRammler
|
||||
)
|
||||
:
|
||||
distributionModel(p),
|
||||
minValue_(p.minValue_),
|
||||
maxValue_(p.maxValue_),
|
||||
d_(p.d_),
|
||||
lambda_(p.lambda_),
|
||||
n_(p.n_)
|
||||
{}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
|
||||
|
||||
Foam::distributionModels::massRosinRammler::~massRosinRammler()
|
||||
{}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
|
||||
|
||||
Foam::scalar Foam::distributionModels::massRosinRammler::sample() const
|
||||
{
|
||||
scalar d;
|
||||
|
||||
// Re-sample if the calculated d is out of the physical range
|
||||
scalar d = 0;
|
||||
do
|
||||
{
|
||||
const scalar a = 3/n_ + 1;
|
||||
const scalar P = rndGen_.sample01<scalar>();
|
||||
const scalar x = Math::invIncGamma(a, P);
|
||||
d = d_*pow(x, 1/n_);
|
||||
} while (d < minValue_ || d > maxValue_);
|
||||
// (YHD:Inverse of Eq. 10)
|
||||
const scalar a = scalar(3)/n_ + scalar(1);
|
||||
const scalar cdfA = Math::incGamma_P(a, pow(minValue_/lambda_, n_) );
|
||||
const scalar cdfB = Math::incGamma_P(a, pow(maxValue_/lambda_, n_) );
|
||||
|
||||
const scalar u = rndGen_.position<scalar>(cdfA, cdfB);
|
||||
const scalar x = Math::invIncGamma(a, u);
|
||||
d = lambda_*pow(x, scalar(1)/n_);
|
||||
|
||||
} while (std::isnan(d));
|
||||
|
||||
return d;
|
||||
}
|
||||
|
||||
|
||||
Foam::scalar Foam::distributionModels::massRosinRammler::minValue() const
|
||||
{
|
||||
return minValue_;
|
||||
}
|
||||
|
||||
|
||||
Foam::scalar Foam::distributionModels::massRosinRammler::maxValue() const
|
||||
{
|
||||
return maxValue_;
|
||||
}
|
||||
|
||||
|
||||
Foam::scalar Foam::distributionModels::massRosinRammler::meanValue() const
|
||||
{
|
||||
return d_;
|
||||
// (YHD:Eqs. 11-12)
|
||||
return lambda_*tgamma(scalar(1)/n_ + scalar(1));
|
||||
}
|
||||
|
||||
|
||||
|
||||
@ -6,6 +6,7 @@
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
Copyright (C) 2016 OpenFOAM Foundation
|
||||
Copyright (C) 2021 OpenCFD Ltd.
|
||||
-------------------------------------------------------------------------------
|
||||
License
|
||||
This file is part of OpenFOAM.
|
||||
@ -24,27 +25,141 @@ License
|
||||
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
|
||||
|
||||
Class
|
||||
Foam::massRosinRammler
|
||||
Foam::distributionModels::massRosinRammler
|
||||
|
||||
Description
|
||||
Mass-based Rosin-Rammler distributionModel.
|
||||
Particle-size distribution model wherein random samples are drawn
|
||||
from the two-parameter Rosin-Rammler (Weibull) probability density
|
||||
function corrected to take into account varying number of particles
|
||||
per parcels for fixed-mass parcels.
|
||||
|
||||
Corrected form of the Rosin-Rammler distribution taking into account the
|
||||
varying number of particels per parces for for fixed-mass parcels. This
|
||||
distribution should be used when
|
||||
"There is a factor of \f$x^3\f$ difference between the size distribution
|
||||
probability density function of individual particles and modeled parcels
|
||||
of particles, \f$ f_{parcels}(D) = x^3 f_{particles}(D) \f$, because the
|
||||
submodel parameter, \f$ PPP \f$ (number of particles per parcel) is based
|
||||
on a fixed mass per parcel which weights the droplet distribution by
|
||||
a factor proportional to \f$ 1/x^3 \f$ (i.e. \f$ PPP = \dot{m}_{total}
|
||||
\Delta_{t_{inj(per-parcel)}} / {m}_{particle} \f$)." (YHD:p. 1374-1375).
|
||||
|
||||
\c massRosinRammler should be preferred over \c RosinRammler
|
||||
when \c parcelBasisType is based on \c mass:
|
||||
\verbatim
|
||||
parcelBasisType mass;
|
||||
\endverbatim
|
||||
|
||||
See equation 10 in reference:
|
||||
The doubly-truncated mass-corrected
|
||||
Rosin-Rammler probability density function:
|
||||
|
||||
\f[
|
||||
f(x; \lambda, n, A, B) =
|
||||
x^3
|
||||
\frac{n}{\lambda}
|
||||
\left( \frac{x}{\lambda} \right)^{n-1}
|
||||
\frac{
|
||||
\exp\{ -(\frac{x}{\lambda})^n \}
|
||||
}
|
||||
{
|
||||
\exp\{ -(\frac{A}{\lambda})^n \}
|
||||
- \exp\{ -(\frac{B}{\lambda})^n \}
|
||||
}
|
||||
\f]
|
||||
where
|
||||
|
||||
\vartable
|
||||
f(x; \lambda, n, A, B) | Rosin-Rammler probability density function
|
||||
\lambda | Scale parameter
|
||||
n | Shape parameter
|
||||
x | Sample
|
||||
A | Minimum of the distribution
|
||||
B | Maximum of the distribution
|
||||
\endvartable
|
||||
|
||||
Constraints:
|
||||
- \f$ \lambda > 0 \f$
|
||||
- \f$ n > 0 \f$
|
||||
|
||||
Random samples are generated by the inverse transform sampling technique
|
||||
by using the quantile function of the doubly-truncated two-parameter
|
||||
mass-corrected Rosin-Rammler (Weibull) probability density function:
|
||||
|
||||
\f[
|
||||
d = \lambda \, Q(a, u)^{1/n}
|
||||
\f]
|
||||
with
|
||||
|
||||
\f[
|
||||
a = \frac{3}{n} + 1
|
||||
\f]
|
||||
where \f$ Q(z_1, z_2) \f$ is the inverse of regularised lower incomplete
|
||||
gamma function, and \f$ u \f$ is a sample drawn from the uniform
|
||||
probability density function on the interval \f$ (x, y) \f$:
|
||||
|
||||
\f[
|
||||
x = \gamma \left( a, \frac{A}{\lambda}^n \right)
|
||||
\f]
|
||||
|
||||
\f[
|
||||
y = \gamma \left( a, \frac{B}{\lambda}^n \right)
|
||||
\f]
|
||||
where \f$ \gamma(z_1, z_2) \f$ is the lower incomplete gamma function.
|
||||
|
||||
Reference:
|
||||
\verbatim
|
||||
Yoon, S. S., Hewson, J. C., DesJardin, P. E., Glaze, D. J.,
|
||||
Black, A. R., & Skaggs, R. R. (2004).
|
||||
Numerical modeling and experimental measurements of a high speed
|
||||
solid-cone water spray for use in fire suppression applications.
|
||||
International Journal of Multiphase Flow, 30(11), 1369-1388.
|
||||
Standard model (tag:YHD):
|
||||
Yoon, S. S., Hewson, J. C., DesJardin, P. E.,
|
||||
Glaze, D. J., Black, A. R., & Skaggs, R. R. (2004).
|
||||
Numerical modeling and experimental measurements of a high speed
|
||||
solid-cone water spray for use in fire suppression applications.
|
||||
International Journal of Multiphase Flow, 30(11), 1369-1388.
|
||||
DOI:10.1016/j.ijmultiphaseflow.2004.07.006
|
||||
\endverbatim
|
||||
|
||||
Usage
|
||||
Minimal example by using \c constant/\<CloudProperties\>:
|
||||
\verbatim
|
||||
subModels
|
||||
{
|
||||
injectionModels
|
||||
{
|
||||
<name>
|
||||
{
|
||||
...
|
||||
parcelBasisType mass;
|
||||
...
|
||||
|
||||
sizeDistribution
|
||||
{
|
||||
type massRosinRammler;
|
||||
massRosinRammlerDistribution
|
||||
{
|
||||
lambda <scaleParameterValue>;
|
||||
n <shapeParameterValue>;
|
||||
minValue <minValue>;
|
||||
maxValue <maxValue>;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
\endverbatim
|
||||
|
||||
where the entries mean:
|
||||
\table
|
||||
Property | Description | Type | Reqd | Deflt
|
||||
type | Type name: massRosinRammler | word | yes | -
|
||||
massRosinRammlerDistribution | Distribution settings | dict | yes | -
|
||||
lambda | Scale parameter | scalar | yes | -
|
||||
n | Shape parameter | scalar | yes | -
|
||||
minValue | Minimum of the distribution | scalar | yes | -
|
||||
maxValue | Maximum of the distribution | scalar | yes | -
|
||||
\endtable
|
||||
|
||||
Note
|
||||
- In the current context, \c lambda (or \c d) is called "characteristic
|
||||
droplet size", and \c n "empirical dimensionless constant to specify
|
||||
the distribution width, sometimes referred to as the dispersion
|
||||
coefficient." (YHD:p. 1374).
|
||||
|
||||
SourceFiles
|
||||
massRosinRammler.C
|
||||
|
||||
@ -70,19 +185,12 @@ class massRosinRammler
|
||||
:
|
||||
public distributionModel
|
||||
{
|
||||
// Private data
|
||||
// Private Data
|
||||
|
||||
//- Distribution minimum
|
||||
scalar minValue_;
|
||||
//- Scale parameter
|
||||
scalar lambda_;
|
||||
|
||||
//- Distribution maximum
|
||||
scalar maxValue_;
|
||||
|
||||
//- Characteristic droplet size
|
||||
scalar d_;
|
||||
|
||||
//- Empirical dimensionless constant to specify the distribution width,
|
||||
// sometimes referred to as the dispersion coefficient
|
||||
//- Shape parameter
|
||||
scalar n_;
|
||||
|
||||
|
||||
@ -97,7 +205,7 @@ public:
|
||||
//- Construct from components
|
||||
massRosinRammler(const dictionary& dict, Random& rndGen);
|
||||
|
||||
//- Construct copy
|
||||
//- Copy construct
|
||||
massRosinRammler(const massRosinRammler& p);
|
||||
|
||||
//- Construct and return a clone
|
||||
@ -106,23 +214,20 @@ public:
|
||||
return autoPtr<distributionModel>(new massRosinRammler(*this));
|
||||
}
|
||||
|
||||
//- No copy assignment
|
||||
void operator=(const massRosinRammler&) = delete;
|
||||
|
||||
|
||||
//- Destructor
|
||||
virtual ~massRosinRammler();
|
||||
virtual ~massRosinRammler() = default;
|
||||
|
||||
|
||||
// Member Functions
|
||||
|
||||
//- Sample the distributionModel
|
||||
//- Sample the distribution
|
||||
virtual scalar sample() const;
|
||||
|
||||
//- Return the minimum value
|
||||
virtual scalar minValue() const;
|
||||
|
||||
//- Return the maximum value
|
||||
virtual scalar maxValue() const;
|
||||
|
||||
//- Return the mean value
|
||||
//- Return the theoretical mean of the distribution
|
||||
virtual scalar meanValue() const;
|
||||
};
|
||||
|
||||
|
||||
@ -6,6 +6,7 @@
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
Copyright (C) 2011-2016 OpenFOAM Foundation
|
||||
Copyright (C) 2021 OpenCFD Ltd.
|
||||
-------------------------------------------------------------------------------
|
||||
License
|
||||
This file is part of OpenFOAM.
|
||||
@ -26,6 +27,9 @@ License
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
#include "multiNormal.H"
|
||||
#include "mathematicalConstants.H"
|
||||
#include "MathFunctions.H"
|
||||
#include "ListOps.H"
|
||||
#include "addToRunTimeSelectionTable.H"
|
||||
|
||||
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
|
||||
@ -48,37 +52,67 @@ Foam::distributionModels::multiNormal::multiNormal
|
||||
)
|
||||
:
|
||||
distributionModel(typeName, dict, rndGen),
|
||||
minValue_(distributionModelDict_.get<scalar>("minValue")),
|
||||
maxValue_(distributionModelDict_.get<scalar>("maxValue")),
|
||||
range_(maxValue_ - minValue_),
|
||||
expectation_(distributionModelDict_.lookup("expectation")),
|
||||
variance_(distributionModelDict_.lookup("variance")),
|
||||
strength_(distributionModelDict_.lookup("strength"))
|
||||
mu_
|
||||
(
|
||||
distributionModelDict_.lookupCompat
|
||||
(
|
||||
"mu",
|
||||
{{"expectation", 2112}}
|
||||
)
|
||||
),
|
||||
sigma_
|
||||
(
|
||||
distributionModelDict_.lookupCompat
|
||||
(
|
||||
"sigma",
|
||||
{{"variance", 2112}}
|
||||
)
|
||||
),
|
||||
weight_
|
||||
(
|
||||
distributionModelDict_.lookupCompat
|
||||
(
|
||||
"weight",
|
||||
{{"strength", 2112}}
|
||||
)
|
||||
)
|
||||
{
|
||||
check();
|
||||
|
||||
scalar sMax = 0;
|
||||
label n = strength_.size();
|
||||
for (label i=0; i<n; i++)
|
||||
scalar sum = 0;
|
||||
for (label i = 0; i < weight_.size(); ++i)
|
||||
{
|
||||
scalar x = expectation_[i];
|
||||
scalar s = strength_[i];
|
||||
for (label j=0; j<n; j++)
|
||||
if (i > 0 && weight_[i] < weight_[i-1])
|
||||
{
|
||||
if (i!=j)
|
||||
{
|
||||
scalar x2 = (x - expectation_[j])/variance_[j];
|
||||
scalar y = exp(-0.5*x2*x2);
|
||||
s += strength_[j]*y;
|
||||
}
|
||||
FatalErrorInFunction
|
||||
<< type() << "distribution: "
|
||||
<< "Weights must be specified in a monotonic order." << nl
|
||||
<< "Please see the row i = " << i << nl
|
||||
<< "weight[i-1] = " << weight_[i-1] << nl
|
||||
<< "weight[i] = " << weight_[i]
|
||||
<< exit(FatalError);
|
||||
}
|
||||
|
||||
sMax = max(sMax, s);
|
||||
sum += weight_[i];
|
||||
}
|
||||
|
||||
for (label i=0; i<n; i++)
|
||||
if (sum < VSMALL)
|
||||
{
|
||||
strength_[i] /= sMax;
|
||||
FatalErrorInFunction
|
||||
<< type() << "distribution: "
|
||||
<< "The sum of weights cannot be zero." << nl
|
||||
<< "weight = " << weight_
|
||||
<< exit(FatalError);
|
||||
}
|
||||
|
||||
for (label i = 1; i < weight_.size(); ++i)
|
||||
{
|
||||
weight_[i] += weight_[i-1];
|
||||
}
|
||||
|
||||
for (auto& w : weight_)
|
||||
{
|
||||
w /= sum;
|
||||
}
|
||||
}
|
||||
|
||||
@ -86,18 +120,9 @@ Foam::distributionModels::multiNormal::multiNormal
|
||||
Foam::distributionModels::multiNormal::multiNormal(const multiNormal& p)
|
||||
:
|
||||
distributionModel(p),
|
||||
minValue_(p.minValue_),
|
||||
maxValue_(p.maxValue_),
|
||||
range_(p.range_),
|
||||
expectation_(p.expectation_),
|
||||
variance_(p.variance_),
|
||||
strength_(p.strength_)
|
||||
{}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
|
||||
|
||||
Foam::distributionModels::multiNormal::~multiNormal()
|
||||
mu_(p.mu_),
|
||||
sigma_(p.sigma_),
|
||||
weight_(p.weight_)
|
||||
{}
|
||||
|
||||
|
||||
@ -105,54 +130,54 @@ Foam::distributionModels::multiNormal::~multiNormal()
|
||||
|
||||
Foam::scalar Foam::distributionModels::multiNormal::sample() const
|
||||
{
|
||||
scalar y = 0;
|
||||
scalar x = 0;
|
||||
label n = expectation_.size();
|
||||
bool success = false;
|
||||
const scalar u = rndGen_.sample01<scalar>();
|
||||
|
||||
while (!success)
|
||||
for (label i = 0; i < weight_.size(); ++i)
|
||||
{
|
||||
x = minValue_ + range_*rndGen_.sample01<scalar>();
|
||||
y = rndGen_.sample01<scalar>();
|
||||
scalar p = 0.0;
|
||||
|
||||
for (label i=0; i<n; i++)
|
||||
if (weight_[i] > u)
|
||||
{
|
||||
scalar nu = expectation_[i];
|
||||
scalar sigma = variance_[i];
|
||||
scalar s = strength_[i];
|
||||
scalar v = (x - nu)/sigma;
|
||||
p += s*exp(-0.5*v*v);
|
||||
}
|
||||
|
||||
if (y<p)
|
||||
{
|
||||
success = true;
|
||||
return sample(mu_[i], sigma_[i]);
|
||||
}
|
||||
}
|
||||
|
||||
return x;
|
||||
const label last = weight_.size() - 1;
|
||||
|
||||
return sample(mu_[last], sigma_[last]);
|
||||
}
|
||||
|
||||
|
||||
Foam::scalar Foam::distributionModels::multiNormal::minValue() const
|
||||
Foam::scalar Foam::distributionModels::multiNormal::sample
|
||||
(
|
||||
const scalar mu,
|
||||
const scalar sigma
|
||||
) const
|
||||
{
|
||||
return minValue_;
|
||||
}
|
||||
const scalar a = (minValue_ - mu)/sigma;
|
||||
const scalar b = (maxValue_ - mu)/sigma;
|
||||
|
||||
const scalar aPhi = 0.5*(1.0 + erf(a/Foam::sqrt(2.0)));
|
||||
const scalar bPhi = 0.5*(1.0 + erf(b/Foam::sqrt(2.0)));
|
||||
|
||||
Foam::scalar Foam::distributionModels::multiNormal::maxValue() const
|
||||
{
|
||||
return maxValue_;
|
||||
const scalar u = rndGen_.sample01<scalar>();
|
||||
const scalar p = u*(bPhi - aPhi) + aPhi;
|
||||
|
||||
// (B:p. 20-24)
|
||||
const scalar x =
|
||||
mu + sigma*Foam::sqrt(2.0)*Math::erfInv(2.0*p - 1.0);
|
||||
|
||||
// Note: numerical approximation of the inverse function yields slight
|
||||
// inaccuracies
|
||||
|
||||
return min(max(x, minValue_), maxValue_);
|
||||
}
|
||||
|
||||
|
||||
Foam::scalar Foam::distributionModels::multiNormal::meanValue() const
|
||||
{
|
||||
scalar mean = 0.0;
|
||||
forAll(strength_, i)
|
||||
scalar mean = 0;
|
||||
forAll(weight_, i)
|
||||
{
|
||||
mean += strength_[i]*expectation_[i];
|
||||
mean += weight_[i]*mu_[i];
|
||||
}
|
||||
|
||||
return mean;
|
||||
|
||||
@ -6,6 +6,7 @@
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
Copyright (C) 2011-2013 OpenFOAM Foundation
|
||||
Copyright (C) 2021 OpenCFD Ltd.
|
||||
-------------------------------------------------------------------------------
|
||||
License
|
||||
This file is part of OpenFOAM.
|
||||
@ -27,12 +28,149 @@ Class
|
||||
Foam::distributionModels::multiNormal
|
||||
|
||||
Description
|
||||
A multiNormal distribution model
|
||||
Particle-size distribution model wherein random samples are drawn
|
||||
from a mixture of a finite set of doubly-truncated univariate normal
|
||||
probability density functions:
|
||||
|
||||
\f[
|
||||
g (\mathbf{x}; \mathbf{\mu}, \mathbf{\sigma}, A, B) =
|
||||
\sum_i w_i f(x_i; \mu_i, \sigma_i, A, B)
|
||||
\f]
|
||||
with for any distribution:
|
||||
|
||||
\f[
|
||||
f(x; \mu, \sigma, A, B) =
|
||||
\frac{1}{\sigma}
|
||||
\frac{
|
||||
\phi \left( \frac{x - \mu}{\sigma} \right)
|
||||
}{
|
||||
\Phi \left( \frac{B - \mu}{\sigma} \right)
|
||||
- \Phi \left( \frac{A - \mu}{\sigma} \right)
|
||||
}
|
||||
\f]
|
||||
where
|
||||
|
||||
\vartable
|
||||
f(x; \mu, \sigma, A, B) | Doubly-truncated univariate normal distribution
|
||||
\mu | Mean of the parent general normal distribution
|
||||
\sigma | Standard deviation of the parent general normal distribution
|
||||
\phi(\cdot) | General normal probability density function
|
||||
\Phi(\cdot) | General normal cumulative distribution function
|
||||
x | Sample
|
||||
A | Minimum of the distribution (the same for each distribution)
|
||||
B | Maximum of the distribution (the same for each distribution)
|
||||
w_i | Weighting factor
|
||||
\endvartable
|
||||
|
||||
Constraints:
|
||||
- \f$ \infty > B > A > 0 \f$
|
||||
- \f$ x \in [B,A] \f$
|
||||
- \f$ \sigma^2 > 0 \f$
|
||||
- \f$ w_i >= 0 \f$
|
||||
|
||||
Random samples are generated by a combination of the inverse transform
|
||||
sampling technique and categorical sampling in three steps:
|
||||
- Draw a sample from the uniform probability density function
|
||||
on the unit interval \f$u = (0, 1)\f$
|
||||
- Find the interval among normalised cumulative weight intervals
|
||||
wherein \f$ u \f$ resides
|
||||
- Draw a sample from the distribution corresponding to the interval by
|
||||
using the quantile function of the doubly-truncated univariate normal
|
||||
probability density function by the following expressions (similar to
|
||||
\c distributionModels::normal):
|
||||
|
||||
\f[
|
||||
x = \mu + \sigma \sqrt{2} \, {erf}^{-1} \left( 2 p - 1 \right)
|
||||
\f]
|
||||
with
|
||||
|
||||
\f[
|
||||
p = u \,
|
||||
\left(
|
||||
\Phi\left(
|
||||
\frac{B - \mu}{\sigma}
|
||||
\right)
|
||||
- \Phi\left(
|
||||
\frac{A - \mu}{\sigma}
|
||||
\right)
|
||||
\right)
|
||||
+ \Phi\left( \frac{A - \mu}{\sigma} \right)
|
||||
\f]
|
||||
|
||||
\f[
|
||||
\Phi(\xi) =
|
||||
\frac{1}{2}
|
||||
\left(
|
||||
1 + {erf}\left(\frac{\xi - \mu}{\sigma \sqrt{2} }\right)
|
||||
\right)
|
||||
\f]
|
||||
where \f$ u \f$ is another sample drawn from the uniform probability
|
||||
density function on the unit interval \f$ (0, 1) \f$.
|
||||
|
||||
Usage
|
||||
Minimal example by using \c constant/\<CloudProperties\>:
|
||||
\verbatim
|
||||
model = sum_i strength_i * exp(-0.5*((x - expectation_i)/variance_i)^2 )
|
||||
subModels
|
||||
{
|
||||
injectionModels
|
||||
{
|
||||
<name>
|
||||
{
|
||||
...
|
||||
|
||||
sizeDistribution
|
||||
{
|
||||
type multiNormal;
|
||||
multiNormalDistribution
|
||||
{
|
||||
minValue <min>;
|
||||
maxValue <max>;
|
||||
mu
|
||||
(
|
||||
<mean1>
|
||||
<mean2>
|
||||
...
|
||||
);
|
||||
sigma
|
||||
(
|
||||
<standard deviation1>
|
||||
<standard deviation2>
|
||||
...
|
||||
);
|
||||
weight
|
||||
(
|
||||
<weight1>
|
||||
<weight2>
|
||||
...
|
||||
);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
\endverbatim
|
||||
|
||||
where the entries mean:
|
||||
\table
|
||||
Property | Description | Type | Reqd | Deflt
|
||||
type | Type name: multiNormal | word | yes | -
|
||||
multiNormalDistribution | Distribution settings | dict | yes | -
|
||||
minValue | Minimum of the distribution | scalar | yes | -
|
||||
maxValue | Maximum of the distribution | scalar | yes | -
|
||||
mu | List of means of parent general normal <!--
|
||||
--> distributions | scalarList | yes | -
|
||||
sigma | List of standard deviations of parent <!--
|
||||
--> general normal distributions | scalarList | yes | -
|
||||
weight | List of weights of a given distribution in <!--
|
||||
--> the distribution mixture | scalarList | yes | -
|
||||
\endtable
|
||||
|
||||
Notes
|
||||
- The sum of normal distributions (i.e. a mixture distribution) should not
|
||||
be confused with the sum of normally-distributed random variables.
|
||||
- \c minValue and \c maxValue are the same for all distributions
|
||||
in the distribution mixture.
|
||||
- \c weight should always be input in a non-decreasing (i.e. monotonic) order.
|
||||
|
||||
SourceFiles
|
||||
multiNormal.C
|
||||
@ -59,23 +197,17 @@ class multiNormal
|
||||
:
|
||||
public distributionModel
|
||||
{
|
||||
// Private data
|
||||
// Private Data
|
||||
|
||||
//- Distribution minimum
|
||||
scalar minValue_;
|
||||
//- List of means of the parent general normal distributions
|
||||
List<scalar> mu_;
|
||||
|
||||
//- Distribution maximum
|
||||
scalar maxValue_;
|
||||
//- List of standard deviations of
|
||||
//- the parent general normal distributions
|
||||
List<scalar> sigma_;
|
||||
|
||||
//- Distribution range
|
||||
scalar range_;
|
||||
|
||||
|
||||
// Model coefficients
|
||||
|
||||
List<scalar> expectation_;
|
||||
List<scalar> variance_;
|
||||
List<scalar> strength_;
|
||||
//- List of weights of a given distribution in the mixture
|
||||
List<scalar> weight_;
|
||||
|
||||
|
||||
public:
|
||||
@ -89,7 +221,7 @@ public:
|
||||
//- Construct from components
|
||||
multiNormal(const dictionary& dict, Random& rndGen);
|
||||
|
||||
//- Construct copy
|
||||
//- Copy construct
|
||||
multiNormal(const multiNormal& p);
|
||||
|
||||
//- Construct and return a clone
|
||||
@ -98,23 +230,23 @@ public:
|
||||
return autoPtr<distributionModel>(new multiNormal(*this));
|
||||
}
|
||||
|
||||
//- No copy assignment
|
||||
void operator=(const multiNormal&) = delete;
|
||||
|
||||
|
||||
//- Destructor
|
||||
virtual ~multiNormal();
|
||||
virtual ~multiNormal() = default;
|
||||
|
||||
|
||||
// Member Functions
|
||||
|
||||
//- Sample the distributionModel
|
||||
//- Sample the distribution
|
||||
virtual scalar sample() const;
|
||||
|
||||
//- Return the minimum value
|
||||
virtual scalar minValue() const;
|
||||
//- Sample the normal distribution
|
||||
scalar sample(const scalar mu, const scalar sigma) const;
|
||||
|
||||
//- Return the maximum value
|
||||
virtual scalar maxValue() const;
|
||||
|
||||
//- Return the mean value
|
||||
//- Return the theoretical mean of the distribution
|
||||
virtual scalar meanValue() const;
|
||||
};
|
||||
|
||||
|
||||
@ -27,8 +27,9 @@ License
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
#include "normal.H"
|
||||
#include "addToRunTimeSelectionTable.H"
|
||||
#include "mathematicalConstants.H"
|
||||
#include "MathFunctions.H"
|
||||
#include "addToRunTimeSelectionTable.H"
|
||||
|
||||
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
|
||||
|
||||
@ -50,44 +51,40 @@ Foam::distributionModels::normal::normal
|
||||
)
|
||||
:
|
||||
distributionModel(typeName, dict, rndGen),
|
||||
minValue_(distributionModelDict_.get<scalar>("minValue")),
|
||||
maxValue_(distributionModelDict_.get<scalar>("maxValue")),
|
||||
expectation_(distributionModelDict_.get<scalar>("expectation")),
|
||||
variance_(distributionModelDict_.get<scalar>("variance")),
|
||||
a_(0.147)
|
||||
mu_
|
||||
(
|
||||
distributionModelDict_.getCompat<scalar>
|
||||
(
|
||||
"mu",
|
||||
{{"expectation", 2112}}
|
||||
)
|
||||
),
|
||||
sigma_
|
||||
(
|
||||
distributionModelDict_.getCompat<scalar>
|
||||
(
|
||||
"sigma",
|
||||
{{"variance", 2112}}
|
||||
)
|
||||
)
|
||||
{
|
||||
if (minValue_ < 0)
|
||||
if (mag(sigma_) == 0)
|
||||
{
|
||||
FatalErrorInFunction
|
||||
<< "Minimum value must be greater than zero. "
|
||||
<< "Supplied minValue = " << minValue_
|
||||
<< abort(FatalError);
|
||||
<< "Standard deviation cannot be zero." << nl
|
||||
<< " sigma = " << sigma_ << nl
|
||||
<< exit(FatalError);
|
||||
}
|
||||
|
||||
if (maxValue_ < minValue_)
|
||||
{
|
||||
FatalErrorInFunction
|
||||
<< "Maximum value is smaller than the minimum value:"
|
||||
<< " maxValue = " << maxValue_ << ", minValue = " << minValue_
|
||||
<< abort(FatalError);
|
||||
}
|
||||
check();
|
||||
}
|
||||
|
||||
|
||||
Foam::distributionModels::normal::normal(const normal& p)
|
||||
:
|
||||
distributionModel(p),
|
||||
minValue_(p.minValue_),
|
||||
maxValue_(p.maxValue_),
|
||||
expectation_(p.expectation_),
|
||||
variance_(p.variance_),
|
||||
a_(p.a_)
|
||||
{}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
|
||||
|
||||
Foam::distributionModels::normal::~normal()
|
||||
mu_(p.mu_),
|
||||
sigma_(p.sigma_)
|
||||
{}
|
||||
|
||||
|
||||
@ -95,37 +92,45 @@ Foam::distributionModels::normal::~normal()
|
||||
|
||||
Foam::scalar Foam::distributionModels::normal::sample() const
|
||||
{
|
||||
const scalar a = (minValue_ - mu_)/sigma_;
|
||||
const scalar b = (maxValue_ - mu_)/sigma_;
|
||||
|
||||
scalar a = erf((minValue_ - expectation_)/variance_);
|
||||
scalar b = erf((maxValue_ - expectation_)/variance_);
|
||||
const scalar aPhi = 0.5*(scalar(1) + erf(a/Foam::sqrt(scalar(2))));
|
||||
const scalar bPhi = 0.5*(scalar(1) + erf(b/Foam::sqrt(scalar(2))));
|
||||
|
||||
scalar y = rndGen_.sample01<scalar>();
|
||||
scalar x = Math::erfInv(y*(b - a) + a)*variance_ + expectation_;
|
||||
const scalar u = rndGen_.sample01<scalar>();
|
||||
const scalar p = u*(bPhi - aPhi) + aPhi;
|
||||
|
||||
// (B:p. 20-24)
|
||||
const scalar x =
|
||||
mu_ + sigma_*Foam::sqrt(scalar(2))*Math::erfInv(scalar(2)*p - scalar(1));
|
||||
|
||||
// Note: numerical approximation of the inverse function yields slight
|
||||
// inaccuracies
|
||||
|
||||
x = min(max(x, minValue_), maxValue_);
|
||||
|
||||
return x;
|
||||
}
|
||||
|
||||
|
||||
Foam::scalar Foam::distributionModels::normal::minValue() const
|
||||
{
|
||||
return minValue_;
|
||||
}
|
||||
|
||||
|
||||
Foam::scalar Foam::distributionModels::normal::maxValue() const
|
||||
{
|
||||
return maxValue_;
|
||||
return min(max(x, minValue_), maxValue_);
|
||||
}
|
||||
|
||||
|
||||
Foam::scalar Foam::distributionModels::normal::meanValue() const
|
||||
{
|
||||
return expectation_;
|
||||
const scalar a = (minValue_ - mu_)/sigma_;
|
||||
const scalar b = (maxValue_ - mu_)/sigma_;
|
||||
|
||||
// (B:p. 2)
|
||||
const scalar aphi =
|
||||
scalar(1)/Foam::sqrt(scalar(2)*constant::mathematical::pi)
|
||||
*exp(-0.5*sqr(a));
|
||||
const scalar bphi =
|
||||
scalar(1)/Foam::sqrt(scalar(2)*constant::mathematical::pi)
|
||||
*exp(-0.5*sqr(b));
|
||||
|
||||
// (B:p. 4)
|
||||
const scalar aPhi = 0.5*(scalar(1) + erf(a/Foam::sqrt(scalar(2))));
|
||||
const scalar bPhi = 0.5*(scalar(1) + erf(b/Foam::sqrt(scalar(2))));
|
||||
|
||||
// (B:p. 25)
|
||||
return mu_ - sigma_*(bphi - aphi)/(bPhi - aPhi + VSMALL);
|
||||
}
|
||||
|
||||
|
||||
|
||||
@ -28,13 +28,124 @@ Class
|
||||
Foam::distributionModels::normal
|
||||
|
||||
Description
|
||||
A normal distribution model
|
||||
Particle-size distribution model wherein random samples are drawn
|
||||
from the doubly-truncated univariate normal probability density function:
|
||||
|
||||
\f[
|
||||
f(x; \mu, \sigma, A, B) =
|
||||
\frac{1}{\sigma}
|
||||
\frac{
|
||||
\phi \left( \frac{x - \mu}{\sigma} \right)
|
||||
}{
|
||||
\Phi \left( \frac{B - \mu}{\sigma} \right)
|
||||
- \Phi \left( \frac{A - \mu}{\sigma} \right)
|
||||
}
|
||||
\f]
|
||||
where
|
||||
|
||||
\vartable
|
||||
f(x; \mu, \sigma, A, B) | Doubly-truncated univariate normal distribution
|
||||
\mu | Mean of the parent general normal distribution
|
||||
\sigma | Standard deviation of the parent general normal distribution
|
||||
\phi(\cdot) | General normal probability density function
|
||||
\Phi(\cdot) | General normal cumulative distribution function
|
||||
x | Sample
|
||||
A | Minimum of the distribution
|
||||
B | Maximum of the distribution
|
||||
\endvartable
|
||||
|
||||
Constraints:
|
||||
- \f$ \infty > B > A > 0 \f$
|
||||
- \f$ x \in [B,A] \f$
|
||||
- \f$ \sigma^2 > 0 \f$
|
||||
|
||||
Random samples are generated by the inverse transform sampling technique
|
||||
by using the quantile function of the doubly-truncated univariate normal
|
||||
probability density function:
|
||||
|
||||
\f[
|
||||
x = \mu + \sigma \sqrt{2} \, {erf}^{-1} \left( 2 p - 1 \right)
|
||||
\f]
|
||||
with
|
||||
|
||||
\f[
|
||||
p = u \,
|
||||
\left(
|
||||
\Phi\left(
|
||||
\frac{B - \mu}{\sigma}
|
||||
\right)
|
||||
- \Phi\left(
|
||||
\frac{A - \mu}{\sigma}
|
||||
\right)
|
||||
\right)
|
||||
+ \Phi\left( \frac{A - \mu}{\sigma} \right)
|
||||
\f]
|
||||
|
||||
\f[
|
||||
\Phi(\xi) =
|
||||
\frac{1}{2}
|
||||
\left(
|
||||
1 + {erf}\left(\frac{\xi - \mu}{\sigma \sqrt{2} }\right)
|
||||
\right)
|
||||
\f]
|
||||
where \f$ u \f$ is sample drawn from the uniform probability
|
||||
density function on the unit interval \f$ (0, 1) \f$.
|
||||
|
||||
Reference:
|
||||
\verbatim
|
||||
model = strength * exp(-0.5*((x - expectation)/variance)^2 )
|
||||
Governing expressions (tag:B):
|
||||
Burkardt, J. (2014).
|
||||
The truncated normal distribution.
|
||||
Department of Scientific Computing Website,
|
||||
Florida State University, 1-35.
|
||||
URL:people.sc.fsu.edu/~jburkardt/presentations/truncated_normal.pdf
|
||||
(Retrieved on: 19 Feb 2021)
|
||||
\endverbatim
|
||||
|
||||
strength only has meaning if there's more than one distribution model
|
||||
Usage
|
||||
Minimal example by using \c constant/\<CloudProperties\>:
|
||||
\verbatim
|
||||
subModels
|
||||
{
|
||||
injectionModels
|
||||
{
|
||||
<name>
|
||||
{
|
||||
...
|
||||
|
||||
sizeDistribution
|
||||
{
|
||||
type normal;
|
||||
normalDistribution
|
||||
{
|
||||
mu <mean>;
|
||||
sigma <stardard deviation>;
|
||||
minValue <min>;
|
||||
maxValue <max>;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
\endverbatim
|
||||
|
||||
where the entries mean:
|
||||
\table
|
||||
Property | Description | Type | Reqd | Deflt
|
||||
type | Type name: normal | word | yes | -
|
||||
normalDistribution | Distribution settings | dict | yes | -
|
||||
mu | Mean of the parent general normal distribution <!--
|
||||
--> | scalar | yes | -
|
||||
sigma | Standard deviation of the parent general normal <!--
|
||||
--> distribution | scalar | yes | -
|
||||
minValue | Minimum of the distribution | scalar | yes | -
|
||||
maxValue | Maximum of the distribution | scalar | yes | -
|
||||
\endtable
|
||||
|
||||
Note
|
||||
- The mean and standard deviation of the parent general normal probability
|
||||
distribution function (i.e. input) are not the same with those of the
|
||||
truncated probability distribution function.
|
||||
|
||||
SourceFiles
|
||||
normal.C
|
||||
@ -63,19 +174,11 @@ class normal
|
||||
{
|
||||
// Private Data
|
||||
|
||||
//- Distribution minimum
|
||||
scalar minValue_;
|
||||
//- Mean of the parent general normal distribution
|
||||
scalar mu_;
|
||||
|
||||
//- Distribution maximum
|
||||
scalar maxValue_;
|
||||
|
||||
|
||||
// Model coefficients
|
||||
|
||||
scalar expectation_;
|
||||
scalar variance_;
|
||||
|
||||
scalar a_;
|
||||
//- Standard deviation of the parent general normal distribution
|
||||
scalar sigma_;
|
||||
|
||||
|
||||
public:
|
||||
@ -89,7 +192,7 @@ public:
|
||||
//- Construct from components
|
||||
normal(const dictionary& dict, Random& rndGen);
|
||||
|
||||
//- Construct copy
|
||||
//- Copy construct
|
||||
normal(const normal& p);
|
||||
|
||||
//- Construct and return a clone
|
||||
@ -98,23 +201,20 @@ public:
|
||||
return autoPtr<distributionModel>(new normal(*this));
|
||||
}
|
||||
|
||||
//- No copy assignment
|
||||
void operator=(const normal&) = delete;
|
||||
|
||||
|
||||
//- Destructor
|
||||
virtual ~normal();
|
||||
virtual ~normal() = default;
|
||||
|
||||
|
||||
// Member Functions
|
||||
|
||||
//- Sample the distributionModel
|
||||
//- Sample the distribution
|
||||
virtual scalar sample() const;
|
||||
|
||||
//- Return the minimum value
|
||||
virtual scalar minValue() const;
|
||||
|
||||
//- Return the maximum value
|
||||
virtual scalar maxValue() const;
|
||||
|
||||
//- Return the mean value
|
||||
//- Return the theoretical mean of the distribution
|
||||
virtual scalar meanValue() const;
|
||||
};
|
||||
|
||||
|
||||
@ -6,6 +6,7 @@
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
Copyright (C) 2011-2016 OpenFOAM Foundation
|
||||
Copyright (C) 2021 OpenCFD Ltd.
|
||||
-------------------------------------------------------------------------------
|
||||
License
|
||||
This file is part of OpenFOAM.
|
||||
@ -47,9 +48,7 @@ Foam::distributionModels::uniform::uniform
|
||||
Random& rndGen
|
||||
)
|
||||
:
|
||||
distributionModel(typeName, dict, rndGen),
|
||||
minValue_(distributionModelDict_.get<scalar>("minValue")),
|
||||
maxValue_(distributionModelDict_.get<scalar>("maxValue"))
|
||||
distributionModel(typeName, dict, rndGen)
|
||||
{
|
||||
check();
|
||||
}
|
||||
@ -57,15 +56,7 @@ Foam::distributionModels::uniform::uniform
|
||||
|
||||
Foam::distributionModels::uniform::uniform(const uniform& p)
|
||||
:
|
||||
distributionModel(p),
|
||||
minValue_(p.minValue_),
|
||||
maxValue_(p.maxValue_)
|
||||
{}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
|
||||
|
||||
Foam::distributionModels::uniform::~uniform()
|
||||
distributionModel(p)
|
||||
{}
|
||||
|
||||
|
||||
@ -77,18 +68,6 @@ Foam::scalar Foam::distributionModels::uniform::sample() const
|
||||
}
|
||||
|
||||
|
||||
Foam::scalar Foam::distributionModels::uniform::minValue() const
|
||||
{
|
||||
return minValue_;
|
||||
}
|
||||
|
||||
|
||||
Foam::scalar Foam::distributionModels::uniform::maxValue() const
|
||||
{
|
||||
return maxValue_;
|
||||
}
|
||||
|
||||
|
||||
Foam::scalar Foam::distributionModels::uniform::meanValue() const
|
||||
{
|
||||
return 0.5*(minValue_ + maxValue_);
|
||||
|
||||
@ -6,6 +6,7 @@
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
Copyright (C) 2011-2016 OpenFOAM Foundation
|
||||
Copyright (C) 2021 OpenCFD Ltd.
|
||||
-------------------------------------------------------------------------------
|
||||
License
|
||||
This file is part of OpenFOAM.
|
||||
@ -27,7 +28,69 @@ Class
|
||||
Foam::distributionModels::uniform
|
||||
|
||||
Description
|
||||
Uniform/equally-weighted distribution model
|
||||
Particle-size distribution model wherein random samples are drawn
|
||||
from the doubly-truncated uniform probability density function:
|
||||
|
||||
\f[
|
||||
f(x; A, B) = \frac{1}{B - A}
|
||||
\f]
|
||||
where
|
||||
|
||||
\vartable
|
||||
f(x; A, B) | Doubly-truncated uniform distribution
|
||||
x | Sample
|
||||
A | Minimum of the distribution
|
||||
B | Maximum of the distribution
|
||||
\endvartable
|
||||
|
||||
Constraints:
|
||||
- \f$ \infty > B > A > 0 \f$
|
||||
- \f$ x \in [B,A] \f$
|
||||
- \f$ \sigma^2 > 0 \f$
|
||||
|
||||
Random samples are generated by the inverse transform sampling technique
|
||||
by using the quantile function of the uniform probability density function:
|
||||
|
||||
\f[
|
||||
x = u \, (B - A) + A
|
||||
\f]
|
||||
|
||||
where \f$ u \f$ is sample drawn from the uniform probability
|
||||
density function on the unit interval \f$ (0, 1) \f$.
|
||||
|
||||
Usage
|
||||
Minimal example by using \c constant/\<CloudProperties\>:
|
||||
\verbatim
|
||||
subModels
|
||||
{
|
||||
injectionModels
|
||||
{
|
||||
<name>
|
||||
{
|
||||
...
|
||||
|
||||
sizeDistribution
|
||||
{
|
||||
type uniform;
|
||||
uniformDistribution
|
||||
{
|
||||
minValue <min>;
|
||||
maxValue <max>;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
\endverbatim
|
||||
|
||||
where the entries mean:
|
||||
\table
|
||||
Property | Description | Type | Reqd | Deflt
|
||||
type | Type name: uniform | word | yes | -
|
||||
uniformDistribution | Distribution settings | dict | yes | -
|
||||
minValue | Minimum of the distribution | scalar | yes | -
|
||||
maxValue | Maximum of the distribution | scalar | yes | -
|
||||
\endtable
|
||||
|
||||
SourceFiles
|
||||
uniform.C
|
||||
@ -54,15 +117,6 @@ class uniform
|
||||
:
|
||||
public distributionModel
|
||||
{
|
||||
// Private data
|
||||
|
||||
//- Distribution minimum
|
||||
scalar minValue_;
|
||||
|
||||
//- Distribution maximum
|
||||
scalar maxValue_;
|
||||
|
||||
|
||||
public:
|
||||
|
||||
//- Runtime type information
|
||||
@ -74,7 +128,7 @@ public:
|
||||
//- Construct from components
|
||||
uniform(const dictionary& dict, Random& rndGen);
|
||||
|
||||
//- Construct copy
|
||||
//- Copy construct
|
||||
uniform(const uniform& p);
|
||||
|
||||
//- Construct and return a clone
|
||||
@ -83,23 +137,20 @@ public:
|
||||
return autoPtr<distributionModel>(new uniform(*this));
|
||||
}
|
||||
|
||||
//- No copy assignment
|
||||
void operator=(const uniform&) = delete;
|
||||
|
||||
|
||||
//- Destructor
|
||||
virtual ~uniform();
|
||||
virtual ~uniform() = default;
|
||||
|
||||
|
||||
// Member Functions
|
||||
|
||||
//- Sample the distributionModel
|
||||
//- Sample the distribution
|
||||
virtual scalar sample() const;
|
||||
|
||||
//- Return the minimum value
|
||||
virtual scalar minValue() const;
|
||||
|
||||
//- Return the maximum value
|
||||
virtual scalar maxValue() const;
|
||||
|
||||
//- Return the mean value
|
||||
//- Return the theoretical mean of the distribution
|
||||
virtual scalar meanValue() const;
|
||||
};
|
||||
|
||||
|
||||
Loading…
Reference in New Issue
Block a user