Spheres/Box contact : spheres ori suddenly set to NaN.

[Kneib François]

Hi,

Please see and launch the attached script.
Some spheres are falling into a box because of gravity. The simulation starts well, but when all particles are almost steady, wait for a few seconds and you will see that some of them who are in contact with the bottom wall will disapear. I spent some time to investigate and that's my preliminary conclusions :
- it is fully reproducible
- particles don't really disappear : their orientation is suddenly turned to Quaternion(Nan,Nan,Nan,2*pi). It's not linked with a high torque on them. All other attributes remains okay, and the other particles still can interact with it. If you manually set the orientation to a correct value after it disappears, it will appear again at the right place.
- you can reproduce this with one sphere (uncomment the last line of the script and comment the previous paragraph) : let it fall then apply to it a little angVel : "O.bodies[-1].state.angVel=(0.,0.,0.001)". It will disappear after a few timesteps.
- it's independent on the timestep value.
- it's dependant on the interpenetration of the contact : if you increase the gravity or decrease the young modulus, you will be able to avoid this behaviour.

Does somebody have an idea ? I think it's a rotation calculation issue, in a very specific geometric case.

[Sebastian Pucilowski (sabamacx)]

Confirm reproducibility with your script.

Yade 0.97.0

[Bruno Chareyre (bruno-chareyre)]

I have an idea... maybe.
Do you see a numerical instability with angVel changing sign at each step before it becomes NaN?

[Kneib François (francois-kneib)]

Not really, even after the sphere disappears, angVel remains OK, but always very small. It also not oscillates at timestep timescale. Below some values of angVel at each time step, before and after body 5 disappears.

O.bodies[5].state.angVel[2] :

9.265986930721519e-161
9.198671437378891e-161
7.230838020503688e-161
3.762575491416308e-161
-5.009685389573714e-162
-3.930775072620818e-161
-6.6976705675296116e-161
-8.314102662474639e-161
-8.495241198140243e-161
-6.9233738378205625e-161

#### Body n°5 disappears ####

-3.877830795599762e-161
2.2185263734549886e-163
3.263461437999207e-161
5.983948968740244e-161
7.704273003639118e-161
8.12129667077455e-161
6.948256554330593e-161

This problem seems to occur in a lot of (not specific) cases, and is much more visible in 2d simulations as a missing spheres let a hole. I'm very surprised that nobody have seen this before ...

Thanks for your replies.

[Kneib François (francois-kneib)]

We are looking at the right direction Bruno :
Spheres that doesn't disappear have a minimum value of angVel of about 1e-11
All spheres that disappear have an angVel value of about 1e-161.

To illustrate that, I wait for all spheres to deposit, and before the first one disappears, I make this (the first three values are the walls) :

for b in O.bodies:
    if(abs(b.state.angVel[2])<1e-20):
        print b.id,b.state.angVel[2]

0 0.0
1 0.0
2 0.0
5 9.26598693072e-161
12 -1.31771647554e-149
48 -2.83070071901e-159
116 -1.20280838825e-134
119 -1.42685193682e-153
192 -3.04156472344e-158
202 -9.57516159925e-161
277 -5.13686301693e-158

Then I wait for spheres to disappear for a while and I do :

for b in O.bodies:
    if(isnan(b.state.ori[0])):
        print b.id

5
12
48
116
119
192
202
277

The bug fix looks like a min threshold value ...

[Bruno Chareyre (bruno-chareyre)]

https://github.com/yade/trunk/blob/master/pkg/dem/NewtonIntegrator.cpp#L236
Division by zero anyone?

[Jérôme Duriez (jduriez)]

Is https://github.com/yade/trunk/blob/master/pkg/dem/NewtonIntegrator.cpp#L235 not enough to avoid division by zero ?

[Bruno Chareyre (bruno-chareyre)]

I'll commit fix in a minute.

[Kneib François (francois-kneib)]

Thank you !

[Bruno Chareyre (bruno-chareyre)]

@Jerome
Unfortunately, your hyperlink now shows the fixed version... (and mine too in #6).
[1] points to the hold one.
And no, the fact that axis!=Vector3r::Zero() does not imply that axis.norm()>0, since MIN_REAL² < MIN_REAL.

[1] https://github.com/yade/trunk/blob/6319de511e8118da4ca120ee1725f34060186a2e/pkg/dem/NewtonIntegrator.cpp#L235