ENH: Updates to thermo/pdfs library

- added Niklas' update to Rosin-Rammler
- new fixedValue type
- code clean-up
This commit is contained in:
andy
2010-03-29 17:32:31 +01:00
parent cb7fd50852
commit b3bf830949
18 changed files with 711 additions and 321 deletions

View File

@ -1,17 +1,12 @@
pdf = pdf
uniform = uniform
normal = normal
general = general
exponential = exponential
RosinRammler = RosinRammler
pdf/pdf.C
pdf/newPdf.C
$(pdf)/pdf.C
$(pdf)/newPdf.C
$(uniform)/uniform.C
$(normal)/normal.C
$(general)/general.C
$(exponential)/exponential.C
$(RosinRammler)/RosinRammler.C
exponential/exponential.C
fixedValue/fixedValue.C
general/general.C
multiNormal/multiNormal.C
normal/normal.C
RosinRammler/RosinRammler.C
uniform/uniform.C
LIB = $(FOAM_LIBBIN)/libpdf

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -31,106 +31,52 @@ License
namespace Foam
{
defineTypeNameAndDebug(RosinRammler, 0);
addToRunTimeSelectionTable(pdf, RosinRammler, dictionary);
namespace pdfs
{
defineTypeNameAndDebug(RosinRammler, 0);
addToRunTimeSelectionTable(pdf, RosinRammler, dictionary);
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::RosinRammler::RosinRammler(const dictionary& dict, Random& rndGen)
Foam::pdfs::RosinRammler::RosinRammler(const dictionary& dict, Random& rndGen)
:
pdf(dict, rndGen),
pdfDict_(dict.subDict(typeName + "PDF")),
pdf(typeName, dict, rndGen),
minValue_(readScalar(pdfDict_.lookup("minValue"))),
maxValue_(readScalar(pdfDict_.lookup("maxValue"))),
d_(pdfDict_.lookup("d")),
n_(pdfDict_.lookup("n")),
ls_(d_),
range_(maxValue_ - minValue_)
d_(readScalar(pdfDict_.lookup("d"))),
n_(readScalar(pdfDict_.lookup("n")))
{
if (minValue_<0)
{
FatalErrorIn
(
"RosinRammler::RosinRammler(const dictionary& dict)"
) << " minValue = " << minValue_ << ", it must be >0." << nl
<< abort(FatalError);
}
if (maxValue_<minValue_)
{
FatalErrorIn
(
"RosinRammler::RosinRammler(const dictionary& dict)"
) << " maxValue is smaller than minValue." << nl
<< abort(FatalError);
}
// find max value so that it can be normalized to 1.0
scalar sMax = 0;
label n = d_.size();
for (label i=0; i<n; i++)
{
scalar s = exp(-1.0);
for (label j=0; j<n; j++)
{
if (i!=j)
{
scalar xx = pow(d_[j]/d_[i], n_[j]);
scalar y = xx*exp(-xx);
s += y;
}
}
sMax = max(sMax, s);
}
for (label i=0; i<n; i++)
{
ls_[i] /= sMax;
}
check();
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::RosinRammler::~RosinRammler()
Foam::pdfs::RosinRammler::~RosinRammler()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::scalar Foam::RosinRammler::sample() const
Foam::scalar Foam::pdfs::RosinRammler::sample() const
{
scalar y = 0;
scalar x = 0;
scalar p = 0.0;
do
{
x = minValue_ + range_*rndGen_.scalar01();
y = rndGen_.scalar01();
p = 0.0;
forAll(d_, i)
{
scalar xx = pow(x/d_[i], n_[i]);
p += ls_[i]*xx*exp(-xx);
}
} while (y>p);
scalar K = 1.0 - exp(-pow((maxValue_ - minValue_)/d_, n_));
scalar y = rndGen_.scalar01();
scalar x = minValue_ + d_*::pow(-log(1.0 - y*K), 1.0/n_);
return x;
}
Foam::scalar Foam::RosinRammler::minValue() const
Foam::scalar Foam::pdfs::RosinRammler::minValue() const
{
return minValue_;
}
Foam::scalar Foam::RosinRammler::maxValue() const
Foam::scalar Foam::pdfs::RosinRammler::maxValue() const
{
return maxValue_;
}

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -29,14 +29,12 @@ Description
Rosin-Rammler pdf
@f[
pdf = ( x/d )^n \exp ( -( x/d )^n )
cumulative pdf = (1.0 - exp( -((x - d0)/d)^n ) / (1.0 - exp( -((d1 - d0)/d)^n )
@f]
SourceFiles
RosinRammlerI.H
RosinRammler.C
RosinRammlerIO.C
\*---------------------------------------------------------------------------*/
@ -49,9 +47,11 @@ SourceFiles
namespace Foam
{
namespace pdfs
{
/*---------------------------------------------------------------------------*\
Class RosinRammler Declaration
Class RosinRammler Declaration
\*---------------------------------------------------------------------------*/
class RosinRammler
@ -60,17 +60,16 @@ class RosinRammler
{
// Private data
dictionary pdfDict_;
//- min and max values of the distribution
//- Distribution minimum
scalar minValue_;
//- Distribution maximum
scalar maxValue_;
List<scalar> d_;
List<scalar> n_;
List<scalar> ls_;
// Model coefficients
scalar range_;
scalar d_;
scalar n_;
public:
@ -108,6 +107,7 @@ public:
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace pdfs
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -24,7 +24,6 @@ License
\*---------------------------------------------------------------------------*/
#include "exponential.H"
#include "addToRunTimeSelectionTable.H"
@ -32,92 +31,49 @@ License
namespace Foam
{
defineTypeNameAndDebug(exponential, 0);
addToRunTimeSelectionTable(pdf, exponential, dictionary);
namespace pdfs
{
defineTypeNameAndDebug(exponential, 0);
addToRunTimeSelectionTable(pdf, exponential, dictionary);
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::exponential::exponential(const dictionary& dict, Random& rndGen)
Foam::pdfs::exponential::exponential(const dictionary& dict, Random& rndGen)
:
pdf(dict, rndGen),
pdfDict_(dict.subDict(typeName + "PDF")),
pdf(typeName, dict, rndGen),
minValue_(readScalar(pdfDict_.lookup("minValue"))),
maxValue_(readScalar(pdfDict_.lookup("maxValue"))),
lambda_(pdfDict_.lookup("lambda")),
ls_(lambda_),
range_(maxValue_ - minValue_)
lambda_(readScalar(pdfDict_.lookup("lambda")))
{
if (minValue_<0)
{
FatalErrorIn
(
"exponential::exponential(const dictionary& dict)"
) << " minValue = " << minValue_ << ", it must be >0." << nl
<< abort(FatalError);
}
scalar sMax = 0;
label n = lambda_.size();
for (label i=0; i<n; i++)
{
scalar s = lambda_[i]*exp(-lambda_[i]*minValue_);
for (label j=0; j<n; j++)
{
if (i!=j)
{
scalar y = lambda_[j]*exp(-lambda_[j]*minValue_);
s += y;
}
}
sMax = max(sMax, s);
}
for (label i=0; i<n; i++)
{
ls_[i] /= sMax;
}
check();
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::exponential::~exponential()
Foam::pdfs::exponential::~exponential()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::scalar Foam::exponential::sample() const
Foam::scalar Foam::pdfs::exponential::sample() const
{
scalar y = 0.0;
scalar x = 0.0;
scalar p = 0.0;
do
{
x = minValue_ + range_*rndGen_.scalar01();
y = rndGen_.scalar01();
p = 0.0;
forAll(lambda_, i)
{
p += ls_[i]*exp(-lambda_[i]*x);
}
} while (p>y);
return x;
scalar y = rndGen_.scalar01();
scalar K = exp(-lambda_*maxValue_) - exp(-lambda_*minValue_);
return -(1.0/lambda_)*log(exp(-lambda_*minValue_) + y*K);
}
Foam::scalar Foam::exponential::minValue() const
Foam::scalar Foam::pdfs::exponential::minValue() const
{
return minValue_;
}
Foam::scalar Foam::exponential::maxValue() const
Foam::scalar Foam::pdfs::exponential::maxValue() const
{
return maxValue_;
}

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -29,9 +29,7 @@ Description
exponential pdf
SourceFiles
exponentialI.H
exponential.C
exponentialIO.C
\*---------------------------------------------------------------------------*/
@ -44,6 +42,8 @@ SourceFiles
namespace Foam
{
namespace pdfs
{
/*---------------------------------------------------------------------------*\
Class exponential Declaration
@ -55,16 +55,17 @@ class exponential
{
// Private data
dictionary pdfDict_;
//- min and max values of the distribution
//- Distribution minimum
scalar minValue_;
//- Distribution maximum
scalar maxValue_;
List<scalar> lambda_;
List<scalar> ls_;
scalar range_;
// Model coefficients
scalar lambda_;
public:
@ -101,14 +102,11 @@ public:
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace pdfs
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
//#include "exponentialI.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,76 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "fixedValue.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
namespace pdfs
{
defineTypeNameAndDebug(fixedValue, 0);
addToRunTimeSelectionTable(pdf, fixedValue, dictionary);
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::pdfs::fixedValue::fixedValue(const dictionary& dict, Random& rndGen)
:
pdf(typeName, dict, rndGen),
value_(readScalar(pdfDict_.lookup("value")))
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::pdfs::fixedValue::~fixedValue()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::scalar Foam::pdfs::fixedValue::fixedValue::sample() const
{
return value_;
}
Foam::scalar Foam::pdfs::fixedValue::fixedValue::minValue() const
{
return -VGREAT;
}
Foam::scalar Foam::pdfs::fixedValue::fixedValue::maxValue() const
{
return VGREAT;
}
// ************************************************************************* //

View File

@ -0,0 +1,103 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2010-2010 OpenCFD Ltd.
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Class
Foam::fixedValue
Description
Returns a fixed value
SourceFiles
fixedValue.C
\*---------------------------------------------------------------------------*/
#ifndef pdfFixedValue_H
#define pdfFixedValue_H
#include "pdf.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace pdfs
{
/*---------------------------------------------------------------------------*\
Class fixedValue Declaration
\*---------------------------------------------------------------------------*/
class fixedValue
:
public pdf
{
// Private data
//- Fixed value
scalar value_;
public:
//- Runtime type information
TypeName("fixedValue");
// Constructors
//- Construct from components
fixedValue
(
const dictionary& dict,
Random& rndGen
);
//- Destructor
virtual ~fixedValue();
// Member Functions
//- Sample the pdf
virtual scalar sample() const;
//- Return the minimum value
virtual scalar minValue() const;
//- Return the maximum value
virtual scalar maxValue() const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace pdfs
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -31,84 +31,109 @@ License
namespace Foam
{
defineTypeNameAndDebug(general, 0);
addToRunTimeSelectionTable(pdf, general, dictionary);
namespace pdfs
{
defineTypeNameAndDebug(general, 0);
addToRunTimeSelectionTable(pdf, general, dictionary);
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::general::general(const dictionary& dict, Random& rndGen)
Foam::pdfs::general::general(const dictionary& dict, Random& rndGen)
:
pdf(dict, rndGen),
pdfDict_(dict.subDict(typeName + "PDF")),
pdf(typeName, dict, rndGen),
xy_(pdfDict_.lookup("distribution")),
nEntries_(xy_.size()),
minValue_(xy_[0][0]),
maxValue_(xy_[nEntries_-1][0]),
range_(maxValue_ - minValue_)
integral_(nEntries_)
{
// normalize the pdf
scalar yMax = 0;
check();
for (label i=0; i<nEntries_; i++)
// normalize the cumulative pdf
integral_[0] = 0.0;
for (label i=1; i<nEntries_; i++)
{
yMax = max(yMax, xy_[i][1]);
scalar k = (xy_[i][1] - xy_[i-1][1])/(xy_[i][0] - xy_[i-1][0]);
scalar d = xy_[i-1][1] - k*xy_[i-1][0];
scalar y1 = xy_[i][0]*(0.5*k*xy_[i][0] + d);
scalar y0 = xy_[i-1][0]*(0.5*k*xy_[i-1][0] + d);
scalar area = y1 - y0;
integral_[i] = area + integral_[i-1];
}
for (label i=0; i<nEntries_; i++)
{
xy_[i][1] /= yMax;
xy_[i][1] /= integral_[nEntries_-1];
integral_[i] /= integral_[nEntries_-1];
}
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::general::~general()
Foam::pdfs::general::~general()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::scalar Foam::general::sample() const
Foam::scalar Foam::pdfs::general::sample() const
{
scalar x = 0.0;
scalar y = 0.0;
scalar p = 0.0;
scalar y = rndGen_.scalar01();
do
// find the interval where y is in the table
label n=1;
while (integral_[n] <= y)
{
x = minValue_ + range_*rndGen_.scalar01();
y = rndGen_.scalar01();
n++;
}
bool intervalFound = false;
label i = -1;
while (!intervalFound)
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];
scalar alpha = y + xy_[n-1][0]*(0.5*k*xy_[n-1][0] + d) - integral_[n-1];
scalar x = 0.0;
// if k is small it's a linear equation, otherwise it's 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);
scalar x1 = -0.5*p + sqrtEr;
scalar x2 = -0.5*p - sqrtEr;
if ((x1 >= xy_[n-1][0]) && (x1 <= xy_[n][0]))
{
i++;
if ((x>xy_[i][0]) && (x<xy_[i+1][0]))
{
intervalFound = true;
}
x = x1;
}
p = xy_[i][1]
+ (x - xy_[i][0])
*(xy_[i+1][1] - xy_[i][1])
/(xy_[i+1][0] - xy_[i][0]);
} while (p>y);
else
{
x = x2;
}
}
else
{
x = alpha/d;
}
return x;
}
Foam::scalar Foam::general::minValue() const
Foam::scalar Foam::pdfs::general::minValue() const
{
return minValue_;
}
Foam::scalar Foam::general::maxValue() const
Foam::scalar Foam::pdfs::general::maxValue() const
{
return maxValue_;
}

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -29,9 +29,7 @@ Description
general pdf
SourceFiles
generalI.H
general.C
generalIO.C
\*---------------------------------------------------------------------------*/
@ -44,6 +42,8 @@ SourceFiles
namespace Foam
{
namespace pdfs
{
/*---------------------------------------------------------------------------*\
Class general Declaration
@ -55,8 +55,6 @@ class general
{
// Private data
dictionary pdfDict_;
typedef VectorSpace<Vector<scalar>, scalar, 2> pair;
List<pair> xy_;
@ -67,7 +65,7 @@ class general
scalar minValue_;
scalar maxValue_;
scalar range_;
List<scalar> integral_;
public:
@ -105,6 +103,7 @@ public:
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace pdfs
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

View File

@ -0,0 +1,133 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "multiNormal.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
namespace pdfs
{
defineTypeNameAndDebug(multiNormal, 0);
addToRunTimeSelectionTable(pdf, multiNormal, dictionary);
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::pdfs::multiNormal::multiNormal(const dictionary& dict, Random& rndGen)
:
pdf(typeName, dict, rndGen),
minValue_(readScalar(pdfDict_.lookup("minValue"))),
maxValue_(readScalar(pdfDict_.lookup("maxValue"))),
range_(maxValue_ - minValue_),
expectation_(pdfDict_.lookup("expectation")),
variance_(pdfDict_.lookup("variance")),
strength_(pdfDict_.lookup("strength"))
{
check();
scalar sMax = 0;
label n = strength_.size();
for (label i=0; i<n; i++)
{
scalar x = expectation_[i];
scalar s = strength_[i];
for (label j=0; j<n; j++)
{
if (i!=j)
{
scalar x2 = (x - expectation_[j])/variance_[j];
scalar y = exp(-0.5*x2*x2);
s += strength_[j]*y;
}
}
sMax = max(sMax, s);
}
for (label i=0; i<n; i++)
{
strength_[i] /= sMax;
}
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::pdfs::multiNormal::~multiNormal()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::scalar Foam::pdfs::multiNormal::sample() const
{
scalar y = 0;
scalar x = 0;
label n = expectation_.size();
bool success = false;
while (!success)
{
x = minValue_ + range_*rndGen_.scalar01();
y = rndGen_.scalar01();
scalar p = 0.0;
for (label i=0; i<n; i++)
{
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 x;
}
Foam::scalar Foam::pdfs::multiNormal::minValue() const
{
return minValue_;
}
Foam::scalar Foam::pdfs::multiNormal::maxValue() const
{
return maxValue_;
}
// ************************************************************************* //

View File

@ -0,0 +1,122 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Class
Foam::multiNormal
Description
A multiNormal pdf
@verbatim
pdf = sum_i strength_i * exp(-0.5*((x - expectation_i)/variance_i)^2 )
@endverbatim
SourceFiles
multiNormal.C
\*---------------------------------------------------------------------------*/
#ifndef multiNormal_H
#define multiNormal_H
#include "pdf.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace pdfs
{
/*---------------------------------------------------------------------------*\
Class multiNormal Declaration
\*---------------------------------------------------------------------------*/
class multiNormal
:
public pdf
{
// Private data
//- Distribution minimum
scalar minValue_;
//- Distribution maximum
scalar maxValue_;
//- Distribution range
scalar range_;
// Model coefficients
List<scalar> expectation_;
List<scalar> variance_;
List<scalar> strength_;
public:
//- Runtime type information
TypeName("multiNormal");
// Constructors
//- Construct from components
multiNormal
(
const dictionary& dict,
Random& rndGen
);
//- Destructor
virtual ~multiNormal();
// Member Functions
//- Sample the pdf
virtual scalar sample() const;
//- Return the minimum value
virtual scalar minValue() const;
//- Return the maximum value
virtual scalar maxValue() const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace pdfs
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -26,98 +26,97 @@ License
#include "normal.H"
#include "addToRunTimeSelectionTable.H"
#include "mathematicalConstants.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(normal, 0);
addToRunTimeSelectionTable(pdf, normal, dictionary);
namespace pdfs
{
defineTypeNameAndDebug(normal, 0);
addToRunTimeSelectionTable(pdf, normal, dictionary);
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::normal::normal(const dictionary& dict, Random& rndGen)
Foam::pdfs::normal::normal(const dictionary& dict, Random& rndGen)
:
pdf(dict, rndGen),
pdfDict_(dict.subDict(typeName + "PDF")),
pdf(typeName, dict, rndGen),
minValue_(readScalar(pdfDict_.lookup("minValue"))),
maxValue_(readScalar(pdfDict_.lookup("maxValue"))),
expectation_(pdfDict_.lookup("expectation")),
variance_(pdfDict_.lookup("variance")),
strength_(pdfDict_.lookup("strength")),
range_(maxValue_ - minValue_)
expectation_(readScalar(pdfDict_.lookup("expectation"))),
variance_(readScalar(pdfDict_.lookup("variance"))),
a_(0.147)
{
scalar sMax = 0;
label n = strength_.size();
for (label i=0; i<n; i++)
if (minValue_ < 0)
{
scalar x = expectation_[i];
scalar s = strength_[i];
for (label j=0; j<n; j++)
{
if (i!=j)
{
scalar x2 = (x-expectation_[j])/variance_[j];
scalar y = exp(-0.5*x2*x2);
s += strength_[j]*y;
}
}
sMax = max(sMax, s);
FatalErrorIn("normal::normal(const dictionary&, Random&)")
<< "Minimum value must be greater than zero. "
<< "Supplied minValue = " << minValue_
<< abort(FatalError);
}
for (label i=0; i<n; i++)
if (maxValue_ < minValue_)
{
strength_[i] /= sMax;
FatalErrorIn("normal::normal(const dictionary&, Random&)")
<< "Maximum value is smaller than the minimum value:"
<< " maxValue = " << maxValue_ << ", minValue = " << minValue_
<< abort(FatalError);
}
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::normal::~normal()
Foam::pdfs::normal::~normal()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::scalar Foam::normal::sample() const
Foam::scalar Foam::pdfs::normal::sample() const
{
scalar x = 0.0;
scalar y = 0.0;
scalar p = 0.0;
do
{
x = minValue_ + range_*rndGen_.scalar01();
y = rndGen_.scalar01();
p = 0.0;
scalar a = erf((minValue_ - expectation_)/variance_);
scalar b = erf((maxValue_ - expectation_)/variance_);
forAll(expectation_, i)
{
scalar nu = expectation_[i];
scalar sigma = variance_[i];
scalar s = strength_[i];
scalar v = (x - nu)/sigma;
p += s*exp(-0.5*v*v);
}
} while (p>y);
scalar y = rndGen_.scalar01();
scalar x = erfInv(y*(b - a) + a)*variance_ + expectation_;
// Note: numerical approximation of the inverse function yields slight
// inaccuracies
x = min(max(x, minValue_), maxValue_);
return x;
}
Foam::scalar Foam::normal::minValue() const
Foam::scalar Foam::pdfs::normal::minValue() const
{
return minValue_;
}
Foam::scalar Foam::normal::maxValue() const
Foam::scalar Foam::pdfs::normal::maxValue() const
{
return maxValue_;
}
Foam::scalar Foam::pdfs::normal::erfInv(const scalar y) const
{
scalar k = 2.0/(constant::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)
{
x *= -1.0;
}
return x;
}
// ************************************************************************* //

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -35,9 +35,7 @@ Description
strength only has meaning if there's more than one pdf
SourceFiles
normalI.H
normal.C
normalIO.C
\*---------------------------------------------------------------------------*/
@ -50,6 +48,8 @@ SourceFiles
namespace Foam
{
namespace pdfs
{
/*---------------------------------------------------------------------------*\
Class normal Declaration
@ -61,17 +61,20 @@ class normal
{
// Private data
dictionary pdfDict_;
//- min and max values of the distribution
//- Distribution minimum
scalar minValue_;
//- Distribution maximum
scalar maxValue_;
List<scalar> expectation_;
List<scalar> variance_;
List<scalar> strength_;
scalar range_;
// Model coefficients
scalar expectation_;
scalar variance_;
scalar a_;
public:
@ -104,11 +107,14 @@ public:
//- Return the maximum value
virtual scalar maxValue() const;
scalar erfInv(const scalar y) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace pdfs
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -28,7 +28,7 @@ License
// * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * * //
Foam::autoPtr<Foam::pdf> Foam::pdf::New
Foam::autoPtr<Foam::pdfs::pdf> Foam::pdfs::pdf::New
(
const dictionary& dict,
Random& rndGen
@ -43,7 +43,7 @@ Foam::autoPtr<Foam::pdf> Foam::pdf::New
if (cstrIter == dictionaryConstructorTablePtr_->end())
{
FatalErrorIn("pdf::New(const dictionary&, Random&)")
FatalErrorIn("pdfs::pdf::New(const dictionary&, Random&)")
<< "unknown pdf type " << pdfType << nl << nl
<< "Valid pdf types are:" << nl
<< dictionaryConstructorTablePtr_->sortedToc()

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -30,22 +30,48 @@ License
namespace Foam
{
defineTypeNameAndDebug(pdf, 0);
defineRunTimeSelectionTable(pdf, dictionary);
namespace pdfs
{
defineTypeNameAndDebug(pdf, 0);
defineRunTimeSelectionTable(pdf, dictionary);
}
}
// * * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * //
void Foam::pdfs::pdf::check() const
{
if (minValue() < 0)
{
FatalErrorIn("pdfs::pdf::check() const")
<< type() << "PDF: Minimum value must be greater than zero." << nl
<< "Supplied minValue = " << minValue()
<< abort(FatalError);
}
if (maxValue() < minValue())
{
FatalErrorIn("pdfs::pdf::check() const")
<< type() << "PDF: Maximum value is smaller than the minimum value:"
<< nl << " maxValue = " << maxValue()
<< ", minValue = " << minValue()
<< abort(FatalError);
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::pdf::pdf(const dictionary& dict, Random& rndGen)
Foam::pdfs::pdf::pdf(const word& name, const dictionary& dict, Random& rndGen)
:
dict_(dict),
pdfDict_(dict.subDict(name + "PDF")),
rndGen_(rndGen)
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::pdf::~pdf()
Foam::pdfs::pdf::~pdf()
{}

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -26,37 +26,27 @@ Class
Foam::pdf
Description
A library of runtime-selectable pdf's.
A library of runtime-selectable PDF's.
Returns a sampled value given the expectation (nu) and variance (sigma^2)
Sample of planned pdfs (beta p. 374-375):
- binomial
- geometric
- Poisson
- hypergeometric
- Pascal
- uniform
Current PDF's include:
- exponential
- fixedValue
- general
- multi-normal
- normal
- Gamma
- chi
- t?
- F?
- beta
- Weibull
- Rayleigh
- Cauchy?
- Rosin-Rammler
- uniform
The pdf is tabulated in equidistant nPoints, in an interval.
These values are integrated to obtain the cumulated pdf,
These values are integrated to obtain the cumulated PDF,
which is then used to change the distribution from unifrom to
the actual pdf.
SourceFiles
pdfI.H
pdf.C
pdfIO.C
newPdf.C
\*---------------------------------------------------------------------------*/
@ -71,6 +61,8 @@ SourceFiles
namespace Foam
{
namespace pdfs
{
/*---------------------------------------------------------------------------*\
Class pdf Declaration
@ -83,11 +75,19 @@ protected:
// Protected data
const dictionary& dict_;
//- Coefficients dictionary
const dictionary pdfDict_;
//- Reference to the randmo number generator
Random& rndGen_;
// Protected member functions
//- Check that the PDF is valid
virtual void check() const;
public:
//-Runtime type information
@ -111,7 +111,7 @@ public:
// Constructors
//- Construct from dictionary
pdf(const dictionary& dict, Random& rndGen);
pdf(const word& name, const dictionary& dict, Random& rndGen);
//- Selector
@ -137,6 +137,7 @@ public:
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace pdfs
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -31,43 +31,47 @@ License
namespace Foam
{
defineTypeNameAndDebug(uniform, 0);
addToRunTimeSelectionTable(pdf, uniform, dictionary);
namespace pdfs
{
defineTypeNameAndDebug(uniform, 0);
addToRunTimeSelectionTable(pdf, uniform, dictionary);
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::uniform::uniform(const dictionary& dict, Random& rndGen)
Foam::pdfs::uniform::uniform(const dictionary& dict, Random& rndGen)
:
pdf(dict, rndGen),
pdfDict_(dict.subDict(typeName + "PDF")),
pdf(typeName, dict, rndGen),
minValue_(readScalar(pdfDict_.lookup("minValue"))),
maxValue_(readScalar(pdfDict_.lookup("maxValue"))),
range_(maxValue_ - minValue_)
{}
{
check();
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::uniform::~uniform()
Foam::pdfs::uniform::~uniform()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::scalar Foam::uniform::sample() const
Foam::scalar Foam::pdfs::uniform::sample() const
{
return (minValue_ + rndGen_.scalar01()*range_);
}
Foam::scalar Foam::uniform::minValue() const
Foam::scalar Foam::pdfs::uniform::minValue() const
{
return minValue_;
}
Foam::scalar Foam::uniform::maxValue() const
Foam::scalar Foam::pdfs::uniform::maxValue() const
{
return maxValue_;
}

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -29,9 +29,7 @@ Description
uniform pdf
SourceFiles
uniformI.H
uniform.C
uniformIO.C
\*---------------------------------------------------------------------------*/
@ -44,6 +42,8 @@ SourceFiles
namespace Foam
{
namespace pdfs
{
/*---------------------------------------------------------------------------*\
Class uniform Declaration
@ -55,13 +55,13 @@ class uniform
{
// Private data
dictionary pdfDict_;
//- min and max values of the distribution
//- Distribution minimum
scalar minValue_;
//- Distribution maximum
scalar maxValue_;
// help-variable to avoid always having to calculate max-min
//- Distribution range
scalar range_;
@ -100,6 +100,7 @@ public:
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace pdfs
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //