Rosin Rammler Distribution Mean Value Calculation

Summary

The meanValue() function of the RosinRammler distribution model results in a floating point exception or nan returns.

Steps to reproduce

Call meanValue() on a RosinRammler distribution model.

Example case

This can easily be tested with following code:

#include "fvCFD.H"                 // include basic openFoam classes
#include <set>
#include "RosinRammler.H"

int main(int argc, char *argv[])   // start main loop
{
    Random   rndGen;
    IStringStream is
    (
        "minValue     10E-06;\n"
        "maxValue     20E-06;\n"
        "lambda    1.699E-05;\n"
        "n            8;\n"
    );
    dictionary dict(is);

    distributionModels::RosinRammler distribution(dict,rndGen);
    Info << "Sample: "<<distribution.sample() << endl;
    Info << distribution.meanValue() << endl;
    return 0;
}

What is the current bug behaviour?

The meanValue() function returns nan.

What is the expected correct behavior?

Should return the theoretical mean of the distribution.

Relevant logs and/or images

Environment information

• OpenFOAM version : v2406
• Operating system : ubuntu
• Hardware info : AMD CPU
• Compiler : g++ 11.4.0 GNU compiler

Hi,

With this value of lambda we get parameter a in the meanValue method (RosinRammler.C:103) equal to 58859.151854031777. Then this value is passed to the tgamma function, which returns inf.

Indeed gamma function of the value is: 1.0371689863 x 10^255183.

added community label

Dear Alexey,

Does the lambda value not represent the mean as it was with the Rosin Rammler distribution model for OpenFOAM v2106? There, the parameter was called d, and the mean function was an inline function to return d.

I am not sure if I understand this correctly, but the cited reference states at the beginning that "...where α and β respectively denote the shape and the scale parameters..." Thus, a is the shape and beta the scale parameter. In the RosinRammler class description, lambda is described as the scale parameter, hence beta in the notation of the reference (DOI:10.2139/ssrn.2690255). Which makes a lot of sense as we divide the left or right truncated value by beta. However, in the implementation of OpenFOAM the shape parameter n, corresponding to a in the reference, is used.

I believe there might be a typo.

See Eq. (5) in the reference. It should be in OpenFOAM's notation:

\beta \Rightarrow \lambda~\text{Scale Parameter}

\alpha \Rightarrow n~\text{Shape Parameter}

\frac{\lambda}{e^{-\left(l/\lambda\right)^n}-e^{-\left(r/\lambda\right)^n}}\left[\gamma\left(1/n+1\right)-\gamma\left(1/n+1\right)\right]

In C++ code it would be (OpenFOAM notation)

    const scalar a = scalar(1)/n_ + scalar(1);
    const scalar qMax = pow(maxValue_/lambda_, n_);
    const scalar qMin = pow(minValue_/lambda_, n_);
    const scalar gMax = Math::incGamma_P(a, qMax);
    const scalar gMin = Math::incGamma_P(a, qMin);
 
    return lambda_/(exp(-qMin) - exp(-qMax))*(gMax - gMin);

Just tested this code and it works

Dear @janga,

You are right. Yet, it means that this functionality was never correctly implemented. Before 02/09/2021 the method simply returned d_ (old name of the scale parameter), after 02/09/2021 the implementation always used scale parameter instead of shape parameter everywhere.

Guess, sample method should be also fixed by swapping lambda_ and n_?

assigned to @kuti

The change seems to have come with #1996 (closed). Will tag @kuti on it, but the first glance looks like lambda_ and n_ are indeed swapped.

Hi @janga , Hi @alexey , many thanks for your inspections; much appreciated.

Let me dig the verification cases for this code change and the code. I will let you know in a week or so.

Hi Jan (@janga), Hi Alexey (@alexey),

Although I had realised the existing issue in the Rosin-Rammler implementation four years ago (link), I can't explain how I did implement the correction in the wrong way. My apologies.

I'm pushing the fix today.

mentioned in commit c6220c16

Please let me know if we can close this ticket, or this needs further work.