general.C

Go to the documentation of this file.

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

\*---------------------------------------------------------------------------*/

#include "general.H"
#include "addToRunTimeSelectionTable.H"

// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

namespace Foam
{
namespace distributionModels
{
    defineTypeNameAndDebug(general, 0);
    addToRunTimeSelectionTable(distributionModel, general, dictionary);
}
}

// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

Foam::distributionModels::general::general
(
    const dictionary& dict,
    Random& rndGen
)
:
    distributionModel(typeName, dict, rndGen),
    xy_(distributionModelDict_.lookup("distribution")),
    nEntries_(xy_.size()),
    minValue_(xy_[0][0]),
    maxValue_(xy_[nEntries_-1][0]),
    meanValue_(0.0),
    integral_(nEntries_),
    cumulative_(distributionModelDict_.lookupOrDefault("cumulative", false))
{
    check();

    // 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;

        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;
        }
    }

    // normalise the distribution
    scalar sumArea = integral_.last();

    for (label i=0; i<nEntries_; i++)
    {
        xy_[i][1] /= sumArea;
        integral_[i] /= sumArea;
    }

    meanValue_ /= sumArea;
    meanValue_ = cumulative_ ? (maxValue_ - meanValue_) : meanValue_;

    info();
}


Foam::distributionModels::general::general(const general& p)
:
    distributionModel(p),
    xy_(p.xy_),
    nEntries_(p.nEntries_),
    minValue_(p.minValue_),
    maxValue_(p.maxValue_),
    meanValue_(p.meanValue_),
    integral_(p.integral_),
    cumulative_(p.cumulative_)
{}


// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //

Foam::distributionModels::general::~general()
{}


// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

Foam::scalar Foam::distributionModels::general::sample() const
{
    scalar y = rndGen_.sample01<scalar>();

    // find the interval where y is in the table
    label n=1;
    while (integral_[n] <= y)
    {
        n++;
    }

    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;

    // if k is small it is a linear equation, otherwise it is of second order
    if (mag(k) > small)
    {
        scalar p = 2.0*d/k;
        scalar q = -2.0*alpha/k;
        scalar sqrtEr = sqrt(0.25*p*p - q);

        scalar x1 = -0.5*p + sqrtEr;
        scalar x2 = -0.5*p - sqrtEr;
        if ((x1 >= xy_[n-1][0]) && (x1 <= xy_[n][0]))
        {
            x = x1;
        }
        else
        {
            x = x2;
        }
    }
    else
    {
        x = alpha/d;
    }

    return x;
}


Foam::scalar Foam::distributionModels::general::minValue() const
{
    return minValue_;
}


Foam::scalar Foam::distributionModels::general::maxValue() const
{
    return maxValue_;
}


Foam::scalar Foam::distributionModels::general::meanValue() const
{
    return meanValue_;
}


// ************************************************************************* //