295 lines
8.1 KiB
C
295 lines
8.1 KiB
C
/*---------------------------------------------------------------------------*\
|
|
========= |
|
|
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
|
\\ / O peration |
|
|
\\ / A nd |
|
|
\\/ M anipulation |
|
|
-------------------------------------------------------------------------------
|
|
Copyright (C) 2011-2016 OpenFOAM Foundation
|
|
-------------------------------------------------------------------------------
|
|
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-Distribution
|
|
|
|
Description
|
|
Test the Distribution class
|
|
|
|
Plot normal distribution test in gnuplot using:
|
|
|
|
\verbatim
|
|
normalDistribution(mean, sigma, x) = \
|
|
sqrt(1.0/(2.0*pi*sigma**2))*exp(-(x - mean)**2.0/(2.0*sigma**2))
|
|
|
|
plot normalDistribution(8.5, 2.5, x), "Distribution_scalar_test_x" w p
|
|
\endverbatim
|
|
|
|
\*---------------------------------------------------------------------------*/
|
|
|
|
#include "vector.H"
|
|
#include "labelVector.H"
|
|
#include "tensor.H"
|
|
#include "Distribution.H"
|
|
#include "Random.H"
|
|
#include "dimensionedTypes.H"
|
|
#include "argList.H"
|
|
#include "PstreamReduceOps.H"
|
|
|
|
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
|
|
|
using namespace Foam;
|
|
|
|
int main(int argc, char *argv[])
|
|
{
|
|
#include "setRootCase.H"
|
|
|
|
Random R(918273);
|
|
|
|
{
|
|
// scalar
|
|
label randomDistributionTestSize = 50000000;
|
|
|
|
Distribution<scalar> dS(scalar(5e-2));
|
|
|
|
Info<< nl << "Distribution<scalar>" << nl
|
|
<< "Sampling "
|
|
<< randomDistributionTestSize
|
|
<< " times from GaussNormal distribution."
|
|
<< endl;
|
|
|
|
for (label i = 0; i < randomDistributionTestSize; i++)
|
|
{
|
|
dS.add(2.5*R.GaussNormal<scalar>() + 8.5);
|
|
}
|
|
|
|
Info<< "Mean " << dS.mean() << nl
|
|
<< "Median " << dS.median()
|
|
<< endl;
|
|
|
|
dS.write("Distribution_scalar_test_1");
|
|
|
|
Distribution<scalar> dS2(scalar(1e-2));
|
|
|
|
Info<< nl << "Distribution<scalar>" << nl
|
|
<< "Sampling "
|
|
<< randomDistributionTestSize
|
|
<< " times from GaussNormal distribution."
|
|
<< endl;
|
|
|
|
for (label i = 0; i < randomDistributionTestSize; i++)
|
|
{
|
|
dS2.add(1.5*R.GaussNormal<scalar>() -6.0);
|
|
}
|
|
|
|
Info<< "Mean " << dS2.mean() << nl
|
|
<< "Median " << dS2.median()
|
|
<< endl;
|
|
|
|
dS2.write("Distribution_scalar_test_2");
|
|
|
|
Info<< nl << "Adding previous two Distribution<scalar>" << endl;
|
|
|
|
dS = dS + dS2;
|
|
|
|
dS.write("Distribution_scalar_test_1+2");
|
|
}
|
|
|
|
if (Pstream::parRun())
|
|
{
|
|
// scalar in parallel
|
|
label randomDistributionTestSize = 100000000;
|
|
|
|
Distribution<scalar> dS(scalar(1e-1));
|
|
|
|
Pout<< "Distribution<scalar>" << nl
|
|
<< "Sampling "
|
|
<< randomDistributionTestSize
|
|
<< " times from uniform distribution."
|
|
<< endl;
|
|
|
|
for (label i = 0; i < randomDistributionTestSize; i++)
|
|
{
|
|
dS.add(R.sample01<scalar>() + 10*Pstream::myProcNo());
|
|
}
|
|
|
|
Pout<< "Mean " << dS.mean() << nl
|
|
<< "Median " << dS.median()
|
|
<< endl;
|
|
|
|
reduce(dS, sumOp<Distribution<scalar>>());
|
|
|
|
if (Pstream::master())
|
|
{
|
|
Info<< "Reducing parallel Distribution<scalar>" << nl
|
|
<< "Mean " << dS.mean() << nl
|
|
<< "Median " << dS.median()
|
|
<< endl;
|
|
|
|
dS.write("Distribution_scalar_test_parallel_reduced");
|
|
}
|
|
}
|
|
|
|
{
|
|
// vector
|
|
Distribution<vector> dV(vector(0.1, 0.05, 0.15));
|
|
|
|
label randomDistributionTestSize = 1000000;
|
|
|
|
Info<< nl << "Distribution<vector>" << nl
|
|
<< "Sampling "
|
|
<< randomDistributionTestSize
|
|
<< " times from uniform and GaussNormal distribution."
|
|
<< endl;
|
|
|
|
for (label i = 0; i < randomDistributionTestSize; i++)
|
|
{
|
|
dV.add(R.sample01<vector>());
|
|
|
|
// Adding separate GaussNormal components with component
|
|
// weights
|
|
|
|
dV.add
|
|
(
|
|
vector
|
|
(
|
|
R.GaussNormal<scalar>()*3.0 + 1.5,
|
|
R.GaussNormal<scalar>()*0.25 + 4.0,
|
|
R.GaussNormal<scalar>()*3.0 - 1.5
|
|
),
|
|
vector(1.0, 2.0, 5.0)
|
|
);
|
|
}
|
|
|
|
Info<< "Mean " << dV.mean() << nl
|
|
<< "Median " << dV.median()
|
|
<< endl;
|
|
|
|
dV.write("Distribution_vector_test");
|
|
}
|
|
|
|
// {
|
|
// // labelVector
|
|
// Distribution<labelVector> dLV(labelVector::one*10);
|
|
|
|
// label randomDistributionTestSize = 2000000;
|
|
|
|
// Info<< nl << "Distribution<labelVector>" << nl
|
|
// << "Sampling "
|
|
// << randomDistributionTestSize
|
|
// << " times from uniform distribution."
|
|
// << endl;
|
|
|
|
// for (label i = 0; i < randomDistributionTestSize; i++)
|
|
// {
|
|
// dLV.add
|
|
// (
|
|
// labelVector
|
|
// (
|
|
// R.integer(-1000, 1000),
|
|
// R.integer(-5000, 5000),
|
|
// R.integer(-2000, 7000)
|
|
// )
|
|
// );
|
|
// }
|
|
|
|
// Info<< "Mean " << dLV.mean() << nl
|
|
// << "Median " << dLV.median()
|
|
// << endl;
|
|
|
|
// dLV.write("Distribution_labelVector_test");
|
|
// }
|
|
|
|
{
|
|
// tensor
|
|
Distribution<tensor> dT(tensor::one*1e-2);
|
|
|
|
label randomDistributionTestSize = 2000000;
|
|
|
|
Info<< nl << "Distribution<tensor>" << nl
|
|
<< "Sampling "
|
|
<< randomDistributionTestSize
|
|
<< " times from uniform distribution."
|
|
<< endl;
|
|
|
|
for (label i = 0; i < randomDistributionTestSize; i++)
|
|
{
|
|
dT.add(R.sample01<tensor>());
|
|
}
|
|
|
|
Info<< "Mean " << dT.mean() << nl
|
|
<< "Median " << dT.median()
|
|
<< endl;
|
|
|
|
dT.write("Distribution_tensor_test");
|
|
}
|
|
|
|
{
|
|
// symmTensor
|
|
Distribution<symmTensor> dSyT(symmTensor::one*1e-2);
|
|
|
|
label randomDistributionTestSize = 2000000;
|
|
|
|
Info<< nl << "Distribution<symmTensor>" << nl
|
|
<< "Sampling "
|
|
<< randomDistributionTestSize
|
|
<< " times from uniform distribution."
|
|
<< endl;
|
|
|
|
for (label i = 0; i < randomDistributionTestSize; i++)
|
|
{
|
|
dSyT.add(R.sample01<symmTensor>());
|
|
}
|
|
|
|
Info<< "Mean " << dSyT.mean() << nl
|
|
<< "Median " << dSyT.median()
|
|
<< endl;
|
|
|
|
dSyT.write("Distribution_symmTensor_test");
|
|
}
|
|
|
|
{
|
|
// sphericalTensor
|
|
Distribution<sphericalTensor> dSpT(sphericalTensor::one*1e-2);
|
|
|
|
label randomDistributionTestSize = 50000000;
|
|
|
|
Info<< nl << "Distribution<sphericalTensor>" << nl
|
|
<< "Sampling "
|
|
<< randomDistributionTestSize
|
|
<< " times from uniform distribution."
|
|
<< endl;
|
|
|
|
for (label i = 0; i < randomDistributionTestSize; i++)
|
|
{
|
|
dSpT.add(R.sample01<sphericalTensor>());
|
|
}
|
|
|
|
Info<< "Mean " << dSpT.mean() << nl
|
|
<< "Median " << dSpT.median()
|
|
<< endl;
|
|
|
|
dSpT.write("Distribution_sphericalTensor_test");
|
|
}
|
|
|
|
Info<< nl << "End" << nl << endl;
|
|
|
|
return 0;
|
|
}
|
|
|
|
|
|
// ************************************************************************* //
|