Merge branch 'master' of github.com-OpenFOAM:OpenFOAM/OpenFOAM-dev

This commit is contained in:
Henry Weller
2019-08-16 15:39:36 +01:00
18 changed files with 360 additions and 33 deletions

View File

@ -0,0 +1,3 @@
Test-Distribution2.C
EXE = $(FOAM_USER_APPBIN)/Test-Distribution2

View File

@ -0,0 +1,5 @@
EXE_INC = \
-I$(LIB_SRC)/lagrangian/distributionModels/lnInclude \
EXE_LIBS = \
-ldistributionModels

View File

@ -0,0 +1,105 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2019 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Application
Test-Distribution2
Description
Test the general distributionModel.
\*---------------------------------------------------------------------------*/
#include "Distribution.H"
#include "Random.H"
#include "dimensionedTypes.H"
#include "argList.H"
#include "distributionModel.H"
#include "IFstream.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
using namespace Foam;
int main(int argc, char *argv[])
{
#include "setRootCase.H"
Random R(918273);
dictionary dict(IFstream("testDict")());
Info << nl << "Testing general distribution" << endl;
Info << nl << "Continuous probability density function:" << endl;
autoPtr<distributionModel> dist1
(
distributionModel::New
(
dict.subDict("densityFunction"),
R
)
);
label randomDistributionTestSize = 50000000;
Distribution<scalar> dS(scalar(1e-6));
Info<< nl
<< "Sampling " << randomDistributionTestSize << " times." << endl;
for (label i = 0; i < randomDistributionTestSize; i++)
{
dS.add(dist1->sample());
}
Info<< "Produced mean " << dS.mean() << endl;
dS.write("densityTest");
dS.clear();
Info << nl << "Discrete probability density function:" << endl;
dist1.clear();
dist1 =
distributionModel::New
(
dict.subDict("cumulativeFunction"),
R
);
Info<< nl
<< "Sampling " << randomDistributionTestSize << " times." << endl;
for (label i = 0; i < randomDistributionTestSize; i++)
{
dS.add(dist1->sample());
}
Info<< "Produced mean " << dS.mean() << endl;
dS.write("cumulativeTest");
Info<< nl << "End" << nl << endl;
return 0;
}
// ************************************************************************* //

View File

@ -0,0 +1,24 @@
#!/bin/sh
gnuplot<<EOF
set terminal postscript eps color enhanced font "Helvetica,15"
set output 'comparison.eps'
set multiplot layout 2, 1
set xlabel 'size (m)'
set ylabel 'PDF (-)'
set title 'Density mode test'
plot 'densityTest_' using 1:2 w l t 'produced',\
'densityTest_expected' w l t 'expected'
set key center
set ylabel 'CDF (-)'
set title 'Cumulative mode test'
plot 'cumulativeTest_cumulative_' using 1:2 w l t 'produced',\
'cumulativeTest_expected' w l t 'expected'
EOF
#------------------------------------------------------------------------------

View File

@ -0,0 +1,17 @@
1.00E-05 0
1.50E-05 0.0002792
2.00E-05 0.0019757
2.50E-05 0.0089476
3.00E-05 0.0267279
3.50E-05 0.0614331
4.00E-05 0.116413
4.50E-05 0.193061
5.00E-05 0.289633
5.50E-05 0.401243
6.00E-05 0.521441
6.50E-05 0.635572
7.00E-05 0.743733
7.50E-05 0.839653
8.00E-05 0.911578
8.50E-05 0.96258
9.00E-05 1

View File

@ -0,0 +1,17 @@
1.00E-05 5.00001
1.50E-05 105.600211
2.00E-05 559.001118
2.50E-05 2183.60437
3.00E-05 4797.6096
3.50E-05 8845.41769
4.00E-05 12777.6256
4.50E-05 17344.2347
5.00E-05 20630.6413
5.50E-05 23251.8465
6.00E-05 24006.048
6.50E-05 20835.0417
7.00E-05 21685.4434
7.50E-05 16003.232
8.00E-05 12266.6245
8.50E-05 7765.41553
9.00E-05 6937.61388

View File

@ -0,0 +1,77 @@
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: dev
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object testDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
densityFunction
{
type general;
generalDistribution
{
distribution
(
(10e-06 0.0025)
(15e-06 0.0528)
(20e-06 0.2795)
(25e-06 1.0918)
(30e-06 2.3988)
(35e-06 4.4227)
(40e-06 6.3888)
(45e-06 8.6721)
(50e-06 10.3153)
(55e-06 11.6259)
(60e-06 12.0030)
(65e-06 10.4175)
(70e-06 10.8427)
(75e-06 8.0016)
(80e-06 6.1333)
(85e-06 3.8827)
(90e-06 3.4688)
);
}
}
cumulativeFunction
{
type general;
generalDistribution
{
cumulative true;
distribution
(
(10e-06 0.0000000)
(15e-06 0.0281384)
(20e-06 0.1972235)
(25e-06 0.8949856)
(30e-06 2.6711166)
(35e-06 6.1421180)
(40e-06 11.643361)
(45e-06 19.306838)
(50e-06 28.968245)
(55e-06 40.132642)
(60e-06 52.155796)
(65e-06 63.564077)
(70e-06 74.381959)
(75e-06 83.97055)
(80e-06 91.162850)
(85e-06 96.259317)
(90e-06 100.00000)
);
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2019 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -52,6 +52,7 @@ Foam::distributionModels::RosinRammler::RosinRammler
n_(readScalar(distributionModelDict_.lookup("n")))
{
check();
info();
}

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2019 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -57,6 +57,13 @@ void Foam::distributionModel::check() const
}
void Foam::distributionModel::info() const
{
Info<< " Distribution min: " << minValue() << " max: " << maxValue()
<< " mean: " << meanValue() << endl;
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::distributionModel::distributionModel

View File

@ -80,6 +80,9 @@ protected:
//- Check that the distribution model is valid
virtual void check() const;
//- Print information about the distribution
void info() const;
public:

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2019 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -51,6 +51,7 @@ Foam::distributionModels::exponential::exponential
lambda_(readScalar(distributionModelDict_.lookup("lambda")))
{
check();
info();
}

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2019 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -47,7 +47,9 @@ Foam::distributionModels::fixedValue::fixedValue
:
distributionModel(typeName, dict, rndGen),
value_(readScalar(distributionModelDict_.lookup("value")))
{}
{
info();
}
Foam::distributionModels::fixedValue::fixedValue(const fixedValue& p)

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2019 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -51,35 +51,89 @@ Foam::distributionModels::general::general
minValue_(xy_[0][0]),
maxValue_(xy_[nEntries_-1][0]),
meanValue_(0.0),
integral_(nEntries_)
integral_(nEntries_),
cumulative_(distributionModelDict_.lookupOrDefault("cumulative", false))
{
check();
// normalize the cumulative distributionModel
// Additional sanity checks
if (cumulative_ && xy_[0][1] != 0)
{
FatalErrorInFunction
<< type() << "distribution: "
<< "The cumulative distribution must start from zero."
<< abort(FatalError);
}
for (label i=0; i<nEntries_; i++)
{
if (i > 0 && xy_[i][0] <= xy_[i-1][0])
{
FatalErrorInFunction
<< type() << "distribution: "
<< "The points must be specified in ascending order."
<< abort(FatalError);
}
if (xy_[i][1] < 0)
{
FatalErrorInFunction
<< type() << "distribution: "
<< "The distribution can't be negative."
<< abort(FatalError);
}
if (i > 0 && cumulative_ && xy_[i][1] < xy_[i-1][1])
{
FatalErrorInFunction
<< type() << "distribution: "
<< "Cumulative distribution must be non-decreasing."
<< abort(FatalError);
}
}
// 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.0;
for (label i=1; i<nEntries_; i++)
{
// Integrating k*x+d
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];
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_);
for (label i=0; i<nEntries_; i++)
{
xy_[i][1] /= sumArea;
integral_[i] /= sumArea;
}
meanValue_ /= sumArea;
meanValue_ = cumulative_ ? (maxValue_ - meanValue_) : meanValue_;
info();
}
@ -116,6 +170,11 @@ Foam::scalar Foam::distributionModels::general::sample() const
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];
if (cumulative_)
{
return (y - d)/k;
}
scalar alpha = y + xy_[n-1][0]*(0.5*k*xy_[n-1][0] + d) - integral_[n-1];
scalar x = 0.0;

View File

@ -25,7 +25,13 @@ Class
Foam::distributionModels::general
Description
general distribution model
A general distribution model where the distribution is specified as
(point, value) pairs. By default the values are assumed to represent
a probability density function, but the model also supports specifying a
cumulative distribution function. In both cases it is assumed that the
function is linear between the specified points.
In both modes of operation the values are automatically normalized.
SourceFiles
general.C
@ -58,18 +64,27 @@ class general
typedef VectorSpace<Vector<scalar>, scalar, 2> pair;
//- List of (x, y=f(x)) pairs
List<pair> xy_;
//- Amount of entries in the xy_ list
label nEntries_;
//- Min and max values of the distribution
//- Distribution minimum
scalar minValue_;
//- Distribution maximum
scalar maxValue_;
//- Distribution mean
scalar meanValue_;
//- Values of cumulative distribution function
List<scalar> integral_;
//- Is the distribution specified as cumulative or as density
Switch cumulative_;
public:

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2016-2018 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2016-2019 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -52,6 +52,7 @@ Foam::distributionModels::massRosinRammler::massRosinRammler
n_(readScalar(distributionModelDict_.lookup("n")))
{
check();
info();
}

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2019 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -78,6 +78,8 @@ Foam::distributionModels::multiNormal::multiNormal
{
strength_[i] /= sMax;
}
info();
}

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2019 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -53,21 +53,8 @@ Foam::distributionModels::normal::normal
variance_(readScalar(distributionModelDict_.lookup("variance"))),
a_(0.147)
{
if (minValue_ < 0)
{
FatalErrorInFunction
<< "Minimum value must be greater than zero. "
<< "Supplied minValue = " << minValue_
<< abort(FatalError);
}
if (maxValue_ < minValue_)
{
FatalErrorInFunction
<< "Maximum value is smaller than the minimum value:"
<< " maxValue = " << maxValue_ << ", minValue = " << minValue_
<< abort(FatalError);
}
check();
info();
}

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2019 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -50,6 +50,7 @@ Foam::distributionModels::uniform::uniform
maxValue_(readScalar(distributionModelDict_.lookup("maxValue")))
{
check();
info();
}