src/thermophysicalModels/pdfs/exponential/exponential.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
   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 the
  13     Free Software Foundation; either version 2 of the License, or (at your
  14     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, write to the Free Software Foundation,
  23     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
  24
  25 \*---------------------------------------------------------------------------*/
  26
  27
  28 #include "exponential.H"
  29 #include "addToRunTimeSelectionTable.H"
  30
  31 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  32
  33 namespace Foam
  34 {
  35     defineTypeNameAndDebug(exponential, 0);
  36     addToRunTimeSelectionTable(pdf, exponential, dictionary);
  37 }
  38
  39 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  40
  41 Foam::exponential::exponential(const dictionary& dict, Random& rndGen)
  42 :
  43     pdf(dict, rndGen),
  44     pdfDict_(dict.subDict(typeName + "PDF")),
  45     minValue_(readScalar(pdfDict_.lookup("minValue"))),
  46     maxValue_(readScalar(pdfDict_.lookup("maxValue"))),
  47     lambda_(pdfDict_.lookup("lambda")),
  48     ls_(lambda_),
  49     range_(maxValue_-minValue_)
  50 {
  51     if (minValue_<0)
  52     {
  53         FatalErrorIn
  54         (
  55             "exponential::exponential(const dictionary& dict)"
  56         ) << " minValue = " << minValue_ << ", it must be >0." << abort(FatalError);
  57     }
  58
  59     scalar sMax = 0;
  60     label n = lambda_.size();
  61     for (label i=0; i<n; i++)
  62     {
  63         scalar s = lambda_[i]*exp(-lambda_[i]*minValue_);
  64         for (label j=0; j<n; j++)
  65         {
  66             if (i!=j)
  67             {
  68                 scalar y = lambda_[j]*exp(-lambda_[j]*minValue_);
  69                 s += y;
  70             }
  71         }
  72
  73         sMax = max(sMax, s);
  74     }
  75
  76     for(label i=0; i<n; i++)
  77     {
  78         ls_[i] /= sMax;
  79     }
  80 }
  81
  82
  83 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
  84
  85 Foam::exponential::~exponential()
  86 {}
  87
  88
  89 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  90
  91 Foam::scalar Foam::exponential::sample() const
  92 {
  93     scalar y = 0;
  94     scalar x = 0;
  95     label n = lambda_.size();
  96     bool success = false;
  97
  98     while (!success)
  99     {
 100         x = minValue_ + range_*rndGen_.scalar01();
 101         y = rndGen_.scalar01();
 102         scalar p = 0.0;
 103
 104         for(label i=0; i<n; i++)
 105         {
 106             p += ls_[i]*exp(-lambda_[i]*x);
 107         }
 108
 109         if (y<p)
 110         {
 111             success = true;
 112         }
 113     }
 114
 115     return x;
 116 }
 117
 118
 119 Foam::scalar Foam::exponential::minValue() const
 120 {
 121     return minValue_;
 122 }
 123
 124
 125 Foam::scalar Foam::exponential::maxValue() const
 126 {
 127     return maxValue_;
 128 }
 129
 130
 131 // ************************************************************************* //