yade-dev team mailing list archive

[Václav Šmilauer <eudoxos@xxxxxxxx>]

> Here https://bugs.launchpad.net/yade/+bug/394687/comments/26
> You said:
> "...Just that everybody knows now that slightly different results are
> OK with openMP (different addition order, hence different rounding; if
> damping!=0, that can make very different results) and it is not a
> bug."
> 
> Could you not explain in a couple of words, why? 
> 
> Can we do something with that?

I thought we've discussed that several times.

In finite precision computations, the order of summation can make the
result different at least sinificant places; try this in python (will be
the same in c++):

(1/10.)+(1/13.)+(1/17.)
(1/17.)+(1/13.)+(1/10.)

you will get

0.23574660633484162
0.23574660633484165

In parallel computation, forces generated by interaction are assigned to
bodies in quasi-random order. Therefore summary force F on body can be
different between single-threaded and multi-threaded, and it can also be
different for different runs of multi-threaded computation depending on
exact thread scheduling by the kernel (inpredictable).

Cundall's damping changes acceleration on body a=F/m to 

ai'=ai(1-d*sgn(Fi*(vi+.5*Δt*ai)), i∈{x,y,z}

where F is force on body, vi is velocity at (t+Δt/2), a is acceleration
at (t), d is dimensionless damping coefficient (say 0.4). The sgn()
function is dicontinuous around zero, therefore there is some (small) ε
that in case |(vi+.5*Δt*ai)|<ε, will make the value of sgn()
arbitrarily ±1, leading non-deterministically to one of

 ai'=1.4*ai
 ai'=0.6*ai

Therefore the rounding "error" at least significant places can lead to
immediate macroscopic effect. 

Note: this is different from integration schemes where little change in
initial conditions makes big difference due to high sensitivity of
continuous equations (e.g. fractals); the equation involved here is
really discontinuous.

I will put it to docs somewhere.

v