Error in 'probability_distribution.h'?

Hi all,

I wanted to add a possibility for the fix template/sphere to work in 2 dimensions. So I needed to add respectively the square expectation values in 'probability_distribution.h'.

I think there is an error for the expectancy value of the uniform distribution that you have entered.
For me it will just be: E[x]= (max^2-min^2)/(2*(max-min)) = (max+min)/2
But in the file you have entered: E[x]=[(max+min)/2]^3

Anton

Forums:

Bug Reports for CFDEM®coupling, LIGGGHTS®, and ParScale

apeskhov | Thu, 01/16/2014 - 16:34

Continue to dig in the source

Continue to dig in the source files to obtain sphere distributions in 2d.

In "fix_template_sphere.cpp", line 185 is certainly wrong. Given that it actually will never be executed due to an error condition on the previous line, it must just be removed.

And another remark concerning the constants you enter in the code. I remarked that they are sometime in a lower precision than the double one. This is the case for E in 'probability_distribution.h', and was also the case for Stokes in 'pair_gran_hooke_history.cpp' (I think the last one is fixed in the third version of LIGGGHTS). While this is perhaps not so critical in engineering applications, a scientist that will try to investigate properties near a critical point will certainly not be happy of this "approximations".

apeskhov | Thu, 01/16/2014 - 19:32

Actually isn't there an error

Actually isn't there an error in the uniform mass distribution function?

The probability distribution function of radius r, for an uniform mass is given by:

P(r)=1/(norm*r^3) where norm=(1/r_min^2-1/r_max^2)/2 .

To obtain this distribution from the uniform one you use the inverse transform sampling. For that we need to compute the quantile function of the previous distribution, i.e. the inverse of the probability of a radius in interval [r_min:R]. After direct computations we find that for a random number rn in the interval [0:1] the radius is given by:

R=sqrt(r_min^2/(1-2*norm*r_min^2*rn))

This is different from the equation you use, which I think is once again wrong.