View Issue Details
ID | Project | Category | View Status | Date Submitted | Last Update |
---|---|---|---|---|---|
0003328 | OpenFOAM | Patch | public | 2019-08-14 11:53 | 2019-08-16 11:44 |
Reporter | tniemi | Assigned To | will | ||
Priority | low | Severity | minor | Reproducibility | N/A |
Status | resolved | Resolution | fixed | ||
Product Version | dev | ||||
Fixed in Version | dev | ||||
Summary | 0003328: Patch: improvements to general distributionModel in the lagrangian library | ||||
Description | I have attached a patch which adds optional keyword to general distributionModel which allows to specify the function as cumulative distribution instead of density function. Quite often data is more readily available in cumulative form, for example from measurements made by sieving particles, and it would be more convenient to be able to directly use that as input. The patch also fixes a bug in the calculation of the mean value of the general distribution. Previously, the mean was calculated as the average (unscaled) probability density which is not correct. Also, some additional sanity checks and documentation are added. The sanity checks enforce that a valid distribution must be specified. I would also propose that the distributionModels could write the min, max and mean values to the log file during construction. The reason for this is that it is quite easy to make a mistake when specifying the distributions (especially general) and output to log would make it easier to double check that everything is in order. The patch adds an info-function and corresponding calls to support this. Additionally, there is a test application which checks the functionality provided by this patch. | ||||
Tags | No tags attached. | ||||
|
patch.diff (22,654 bytes)
diff --git a/applications/test/Distribution2/Make/files b/applications/test/Distribution2/Make/files new file mode 100644 index 0000000..723a0b5 --- /dev/null +++ b/applications/test/Distribution2/Make/files @@ -0,0 +1,3 @@ +Test-Distribution2.C + +EXE = $(FOAM_USER_APPBIN)/Test-Distribution2 diff --git a/applications/test/Distribution2/Make/options b/applications/test/Distribution2/Make/options new file mode 100644 index 0000000..ae86ebe --- /dev/null +++ b/applications/test/Distribution2/Make/options @@ -0,0 +1,5 @@ +EXE_INC = \ + -I$(LIB_SRC)/lagrangian/distributionModels/lnInclude \ + +EXE_LIBS = \ + -ldistributionModels diff --git a/applications/test/Distribution2/Test-Distribution2.C b/applications/test/Distribution2/Test-Distribution2.C new file mode 100644 index 0000000..96816a4 --- /dev/null +++ b/applications/test/Distribution2/Test-Distribution2.C @@ -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; +} + + +// ************************************************************************* // diff --git a/applications/test/Distribution2/createGraphs b/applications/test/Distribution2/createGraphs new file mode 100755 index 0000000..bbcb522 --- /dev/null +++ b/applications/test/Distribution2/createGraphs @@ -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 + + +#------------------------------------------------------------------------------ diff --git a/applications/test/Distribution2/cumulativeTest_expected b/applications/test/Distribution2/cumulativeTest_expected new file mode 100644 index 0000000..0974a41 --- /dev/null +++ b/applications/test/Distribution2/cumulativeTest_expected @@ -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 diff --git a/applications/test/Distribution2/densityTest_expected b/applications/test/Distribution2/densityTest_expected new file mode 100644 index 0000000..87959f8 --- /dev/null +++ b/applications/test/Distribution2/densityTest_expected @@ -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 diff --git a/applications/test/Distribution2/testDict b/applications/test/Distribution2/testDict new file mode 100644 index 0000000..8cb3c84 --- /dev/null +++ b/applications/test/Distribution2/testDict @@ -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) + ); + } +} + + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/lagrangian/distributionModels/RosinRammler/RosinRammler.C b/src/lagrangian/distributionModels/RosinRammler/RosinRammler.C index b739076..9d7e37d 100644 --- a/src/lagrangian/distributionModels/RosinRammler/RosinRammler.C +++ b/src/lagrangian/distributionModels/RosinRammler/RosinRammler.C @@ -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(); } diff --git a/src/lagrangian/distributionModels/distributionModel/distributionModel.C b/src/lagrangian/distributionModels/distributionModel/distributionModel.C index bba012b..6367725 100644 --- a/src/lagrangian/distributionModels/distributionModel/distributionModel.C +++ b/src/lagrangian/distributionModels/distributionModel/distributionModel.C @@ -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 diff --git a/src/lagrangian/distributionModels/distributionModel/distributionModel.H b/src/lagrangian/distributionModels/distributionModel/distributionModel.H index 0b830b1..e7ff181 100644 --- a/src/lagrangian/distributionModels/distributionModel/distributionModel.H +++ b/src/lagrangian/distributionModels/distributionModel/distributionModel.H @@ -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: diff --git a/src/lagrangian/distributionModels/exponential/exponential.C b/src/lagrangian/distributionModels/exponential/exponential.C index fd871ee..6db521e 100644 --- a/src/lagrangian/distributionModels/exponential/exponential.C +++ b/src/lagrangian/distributionModels/exponential/exponential.C @@ -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(); } diff --git a/src/lagrangian/distributionModels/fixedValue/fixedValue.C b/src/lagrangian/distributionModels/fixedValue/fixedValue.C index 0d2fef5..79fa7b4 100644 --- a/src/lagrangian/distributionModels/fixedValue/fixedValue.C +++ b/src/lagrangian/distributionModels/fixedValue/fixedValue.C @@ -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) diff --git a/src/lagrangian/distributionModels/general/general.C b/src/lagrangian/distributionModels/general/general.C index 248eab4..c3911dd 100644 --- a/src/lagrangian/distributionModels/general/general.C +++ b/src/lagrangian/distributionModels/general/general.C @@ -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; diff --git a/src/lagrangian/distributionModels/general/general.H b/src/lagrangian/distributionModels/general/general.H index 52ac38c..05c2882 100644 --- a/src/lagrangian/distributionModels/general/general.H +++ b/src/lagrangian/distributionModels/general/general.H @@ -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: diff --git a/src/lagrangian/distributionModels/massRosinRammler/massRosinRammler.C b/src/lagrangian/distributionModels/massRosinRammler/massRosinRammler.C index f79856f..bfac21f 100644 --- a/src/lagrangian/distributionModels/massRosinRammler/massRosinRammler.C +++ b/src/lagrangian/distributionModels/massRosinRammler/massRosinRammler.C @@ -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(); } diff --git a/src/lagrangian/distributionModels/multiNormal/multiNormal.C b/src/lagrangian/distributionModels/multiNormal/multiNormal.C index 5f99be8..e098f84 100644 --- a/src/lagrangian/distributionModels/multiNormal/multiNormal.C +++ b/src/lagrangian/distributionModels/multiNormal/multiNormal.C @@ -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(); } diff --git a/src/lagrangian/distributionModels/normal/normal.C b/src/lagrangian/distributionModels/normal/normal.C index 7c9d3d6..037a87c 100644 --- a/src/lagrangian/distributionModels/normal/normal.C +++ b/src/lagrangian/distributionModels/normal/normal.C @@ -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(); } diff --git a/src/lagrangian/distributionModels/uniform/uniform.C b/src/lagrangian/distributionModels/uniform/uniform.C index 56cc2e6..0aa726b 100644 --- a/src/lagrangian/distributionModels/uniform/uniform.C +++ b/src/lagrangian/distributionModels/uniform/uniform.C @@ -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(); } |
|
Thanks for the patch, Timo. Pushed as commit https://github.com/OpenFOAM/OpenFOAM-dev/commit/1a54b2ecfc6b1503d10690244bb23a3f5ed985bc |
Date Modified | Username | Field | Change |
---|---|---|---|
2019-08-14 11:53 | tniemi | New Issue | |
2019-08-14 11:53 | tniemi | File Added: patch.diff | |
2019-08-16 11:44 | will | Assigned To | => will |
2019-08-16 11:44 | will | Status | new => resolved |
2019-08-16 11:44 | will | Resolution | open => fixed |
2019-08-16 11:44 | will | Fixed in Version | => dev |
2019-08-16 11:44 | will | Note Added: 0010676 |