applications/test/Distribution/Test-Distribution.C

   1 /*---------------------------------------------------------------------------*\
   2  =========                   |
   3  \\      /   F ield          | OpenFOAM: The Open Source CFD Toolbox
   4   \\    /    O peration      |
   5    \\  /     A nd            | Copyright (C) 2011 OpenFOAM Foundation
   6     \\/      M anipulation   |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of OpenFOAM.
  10
  11     OpenFOAM is free software: you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by
  13     the Free Software Foundation, either version 3 of the License, or
  14     (at your option) any later version.
  15
  16     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
  17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  19     for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 Application
  25     Test-Distribution
  26
  27 Description
  28     Test the Distribution class
  29
  30     Plot normal distribution test in gnuplot using:
  31
  32     \verbatim
  33     normalDistribution(mean, sigma, x) = \
  34         sqrt(1.0/(2.0*pi*sigma**2))*exp(-(x - mean)**2.0/(2.0*sigma**2))
  35
  36     plot normalDistribution(8.5, 2.5, x), "Distribution_scalar_test_x" w p
  37     \endverbatim
  38
  39 \*---------------------------------------------------------------------------*/
  40
  41 #include "vector.H"
  42 #include "labelVector.H"
  43 #include "tensor.H"
  44 #include "Distribution.H"
  45 #include "Random.H"
  46 #include "dimensionedTypes.H"
  47 #include "argList.H"
  48 #include "PstreamReduceOps.H"
  49
  50 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  51
  52 using namespace Foam;
  53
  54 int main(int argc, char *argv[])
  55 {
  56     #   include "setRootCase.H"
  57
  58     Random R(918273);
  59
  60     {
  61         // scalar
  62         label randomDistributionTestSize = 50000000;
  63
  64         Distribution<scalar> dS(scalar(5e-2));
  65
  66         Info<< nl << "Distribution<scalar>" << nl
  67             << "Sampling "
  68             << randomDistributionTestSize
  69             << " times from GaussNormal distribution."
  70             << endl;
  71
  72         for (label i = 0; i < randomDistributionTestSize; i++)
  73         {
  74             dS.add(2.5*R.GaussNormal() + 8.5);
  75         }
  76
  77         Info<< "Mean " << dS.mean() << nl
  78             << "Median " << dS.median()
  79             << endl;
  80
  81         dS.write("Distribution_scalar_test_1");
  82
  83         Distribution<scalar> dS2(scalar(1e-2));
  84
  85         Info<< nl << "Distribution<scalar>" << nl
  86             << "Sampling "
  87             << randomDistributionTestSize
  88             << " times from GaussNormal distribution."
  89             << endl;
  90
  91         for (label i = 0; i < randomDistributionTestSize; i++)
  92         {
  93             dS2.add(1.5*R.GaussNormal() -6.0);
  94         }
  95
  96         Info<< "Mean " << dS2.mean() << nl
  97             << "Median " << dS2.median()
  98             << endl;
  99
 100         dS2.write("Distribution_scalar_test_2");
 101
 102         Info<< nl << "Adding previous two Distribution<scalar>" << endl;
 103
 104         dS = dS + dS2;
 105
 106         dS.write("Distribution_scalar_test_1+2");
 107     }
 108
 109     if (Pstream::parRun())
 110     {
 111         // scalar in parallel
 112         label randomDistributionTestSize = 100000000;
 113
 114         Distribution<scalar> dS(scalar(1e-1));
 115
 116         Pout<< "Distribution<scalar>" << nl
 117             << "Sampling "
 118             << randomDistributionTestSize
 119             << " times from uniform distribution."
 120             << endl;
 121
 122         for (label i = 0; i < randomDistributionTestSize; i++)
 123         {
 124             dS.add(R.scalar01() + 10*Pstream::myProcNo());
 125         }
 126
 127         Pout<< "Mean " << dS.mean() << nl
 128             << "Median " << dS.median()
 129             << endl;
 130
 131         reduce(dS, sumOp< Distribution<scalar> >());
 132
 133         if (Pstream::master())
 134         {
 135             Info<< "Reducing parallel Distribution<scalar>" << nl
 136                 << "Mean " << dS.mean() << nl
 137                 << "Median " << dS.median()
 138                 << endl;
 139
 140             dS.write("Distribution_scalar_test_parallel_reduced");
 141         }
 142     }
 143
 144     {
 145         // vector
 146         Distribution<vector> dV(vector(0.1, 0.05, 0.15));
 147
 148         label randomDistributionTestSize = 1000000;
 149
 150         Info<< nl << "Distribution<vector>" << nl
 151             << "Sampling "
 152             << randomDistributionTestSize
 153             << " times from uniform and GaussNormal distribution."
 154             << endl;
 155
 156         for (label i = 0; i < randomDistributionTestSize; i++)
 157         {
 158             dV.add(R.vector01());
 159
 160             // Adding separate GaussNormal components with component
 161             // weights
 162
 163             dV.add
 164             (
 165                 vector
 166                 (
 167                     R.GaussNormal()*3.0 + 1.5,
 168                     R.GaussNormal()*0.25 + 4.0,
 169                     R.GaussNormal()*3.0 - 1.5
 170                 ),
 171                 vector(1.0, 2.0, 5.0)
 172             );
 173         }
 174
 175         Info<< "Mean " << dV.mean() << nl
 176             << "Median " << dV.median()
 177             << endl;
 178
 179         dV.write("Distribution_vector_test");
 180     }
 181
 182     // {
 183     //     // labelVector
 184     //     Distribution<labelVector> dLV(labelVector::one*10);
 185
 186     //     label randomDistributionTestSize = 2000000;
 187
 188     //     Info<< nl << "Distribution<labelVector>" << nl
 189     //         << "Sampling "
 190     //         << randomDistributionTestSize
 191     //         << " times from uniform distribution."
 192     //         << endl;
 193
 194     //     for (label i = 0; i < randomDistributionTestSize; i++)
 195     //     {
 196     //         dLV.add
 197     //         (
 198     //             labelVector
 199     //             (
 200     //                 R.integer(-1000, 1000),
 201     //                 R.integer(-5000, 5000),
 202     //                 R.integer(-2000, 7000)
 203     //             )
 204     //         );
 205     //     }
 206
 207     //     Info<< "Mean " << dLV.mean() << nl
 208     //         << "Median " << dLV.median()
 209     //         << endl;
 210
 211     //     dLV.write("Distribution_labelVector_test");
 212     // }
 213
 214     {
 215         // tensor
 216         Distribution<tensor> dT(tensor::one*1e-2);
 217
 218         label randomDistributionTestSize = 2000000;
 219
 220         Info<< nl << "Distribution<tensor>" << nl
 221             << "Sampling "
 222             << randomDistributionTestSize
 223             << " times from uniform distribution."
 224             << endl;
 225
 226         for (label i = 0; i < randomDistributionTestSize; i++)
 227         {
 228             dT.add(R.tensor01());
 229         }
 230
 231         Info<< "Mean " << dT.mean() << nl
 232             << "Median " << dT.median()
 233             << endl;
 234
 235         dT.write("Distribution_tensor_test");
 236     }
 237
 238     {
 239         // symmTensor
 240         Distribution<symmTensor> dSyT(symmTensor::one*1e-2);
 241
 242         label randomDistributionTestSize = 2000000;
 243
 244         Info<< nl << "Distribution<symmTensor>" << nl
 245             << "Sampling "
 246             << randomDistributionTestSize
 247             << " times from uniform distribution."
 248             << endl;
 249
 250         for (label i = 0; i < randomDistributionTestSize; i++)
 251         {
 252             dSyT.add(R.symmTensor01());
 253         }
 254
 255         Info<< "Mean " << dSyT.mean() << nl
 256             << "Median " << dSyT.median()
 257             << endl;
 258
 259         dSyT.write("Distribution_symmTensor_test");
 260     }
 261
 262     {
 263         // sphericalTensor
 264         Distribution<sphericalTensor> dSpT(sphericalTensor::one*1e-2);
 265
 266         label randomDistributionTestSize = 50000000;
 267
 268         Info<< nl << "Distribution<sphericalTensor>" << nl
 269             << "Sampling "
 270             << randomDistributionTestSize
 271             << " times from uniform distribution."
 272             << endl;
 273
 274         for (label i = 0; i < randomDistributionTestSize; i++)
 275         {
 276             dSpT.add(R.sphericalTensor01());
 277         }
 278
 279         Info<< "Mean " << dSpT.mean() << nl
 280             << "Median " << dSpT.median()
 281             << endl;
 282
 283         dSpT.write("Distribution_sphericalTensor_test");
 284     }
 285
 286     Info<< nl << "End" << nl << endl;
 287
 288     return 0;
 289 }
 290
 291
 292 // ************************************************************************* //