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

lagrangian: distributionModels: Cumulative general distribution

Browse Source 
The general distribution has been extended to accept cumulative
distribution data, by means of a "cumulative" switch. The calculation of
the mean value has also been corrected for this distribution, and
additional header documentation and parameter checking has been added.

In addition, the distribution models now all print some basic
information (min, max and mean) into the log file to help in checking
that the specification is correct.

Patch contributed by Timo Niemi, VTT.
This commit is contained in:
Will Bainbridge
2019-08-16 09:01:30 +01:00
parent 91e00b40b4
commit 1a54b2ecfc
18 changed files with 360 additions and 33 deletions
Download Patch File Download Diff File Expand all files Collapse all files

3
applications/test/Distribution2/Make/files Normal file
View File

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

5
applications/test/Distribution2/Make/options Normal file
View File

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

105
applications/test/Distribution2/Test-Distribution2.C Normal file
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;
}
// ************************************************************************* //

24
applications/test/Distribution2/createGraphs Executable file
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
#------------------------------------------------------------------------------

17
applications/test/Distribution2/cumulativeTest_expected Normal file
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

17
applications/test/Distribution2/densityTest_expected Normal file
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

77
applications/test/Distribution2/testDict Normal file
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)
);
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

3
src/lagrangian/distributionModels/RosinRammler/RosinRammler.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-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();
}

9
src/lagrangian/distributionModels/distributionModel/distributionModel.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-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

3
src/lagrangian/distributionModels/distributionModel/distributionModel.H
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:

3
src/lagrangian/distributionModels/exponential/exponential.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-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();
}

6
src/lagrangian/distributionModels/fixedValue/fixedValue.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-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)

73
src/lagrangian/distributionModels/general/general.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-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;

19
src/lagrangian/distributionModels/general/general.H
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:

3
src/lagrangian/distributionModels/massRosinRammler/massRosinRammler.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 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();
}

4
src/lagrangian/distributionModels/multiNormal/multiNormal.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-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();
}

19
src/lagrangian/distributionModels/normal/normal.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-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();
}

3
src/lagrangian/distributionModels/uniform/uniform.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-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();
}