bug in neverErase?!

[Christian Jakob]

If I set neverErase=True Yade crashes with following output:

O.run(3000,True)
RuntimeError: std::exception

and also:

Traceback (most recent call last):
  File "/home/me/YADE/YADE2970/lib/yade-bzr2970/py/yade/qt/__init__.py", line 192, in refreshEvent 15:01:13
    self.refreshValues() 15:01:30
  File "/home/me/YADE/YADE2970/lib/yade-bzr2970/py/yade/qt/__init__.py", line 265, in refreshValues 15:01:34
    s=int(t); ms=int(t*1000)%1000; us=int(t*1000000)%1000; ns=int(t*1000000000)%1000
ValueError: cannot convert float NaN to integer

You can reproduce the bug with the attached script (similar to the one i submitted with version 2960).
First try to set Law2_ScGeom_MindlinPhys_Mindlin(neverErase=False), this will work.
Then set neverErase=True and yade crashes.

[Anton Gladky (gladky-anton)]

Hi Christian,

thanks for filing the bug-report, the issue needs to be checked.
But you have in your script something like that:

O.bodies[ii].shape.radius=O.bodies[ii].shape.radius*rad_mul[0]

It is potentially dangerous, as was said before.
Could you, please, check the same script without changing the size of
sphere in runtime?

Thanks.
Anton

[Bruno Chareyre (bruno-chareyre)]

Hi,
I suspect the time-step determination fails because kn and ks are not defined for distant interactions when using Hertz functors (with neverErase=false, there are no distant interactions).
Hertz Ip2 should at least assign kn=ks=0 in that case.
Does it work if you remove GSTimestepper and assign fixed dt?

Actually, another fix would be to set default kn and ks =0 instead of the current NaNs...

[Chiara Modenese (chiara-modenese)]

On 28 November 2011 15:13, Chareyre <email address hidden> wrote:

> Hi,
> I suspect the time-step determination fails because kn and ks are not
> defined for distant interactions when using Hertz functors (with
> neverErase=false, there are no distant interactions).
>
Correct.

> Hertz Ip2 should at least assign kn=ks=0 in that case.
> Does it work if you remove GSTimestepper and assign fixed dt?
>
I agree that Christian should use a fixed time. You can estimate it using
PWave time step (
https://www.yade-dem.org/doc/yade.utils.html?highlight=pwave#yade.utils.PWaveTimeStep
).

Chiara

>
> Actually, another fix would be to set default kn and ks =0 instead of
> the current NaNs...
>
> --
> You received this bug notification because you are a member of Yade
> developers, which is the registrant for Yade.
> https://bugs.launchpad.net/bugs/897237
>
> Title:
> bug in neverErase?!
>
> Status in Yet Another Dynamic Engine:
> New
>
> Bug description:
> If I set neverErase=True Yade crashes with following output:
>
> O.run(3000,True)
> RuntimeError: std::exception
>
> and also:
>
> Traceback (most recent call last):
> File "/home/me/YADE/YADE2970/lib/yade-bzr2970/py/yade/qt/__init__.py",
> line 192, in refreshEvent
> 15:01:13
> self.refreshValues()
>
> 15:01:30
> File "/home/me/YADE/YADE2970/lib/yade-bzr2970/py/yade/qt/__init__.py",
> line 265, in refreshValues
> 15:01:34
> s=int(t); ms=int(t*1000)%1000; us=int(t*1000000)%1000;
> ns=int(t*1000000000)%1000
> ValueError: cannot convert float NaN to integer
>
>
> You can reproduce the bug with the attached script (similar to the one i
> submitted with version 2960).
> First try to set Law2_ScGeom_MindlinPhys_Mindlin(neverErase=False), this
> will work.
> Then set neverErase=True and yade crashes.
>
> To manage notifications about this bug go to:
> https://bugs.launchpad.net/yade/+bug/897237/+subscriptions
>
> _______________________________________________
> Mailing list: https://launchpad.net/~yade-dev
> Post to : <email address hidden>
> Unsubscribe : https://launchpad.net/~yade-dev
> More help : https://help.launchpad.net/ListHelp
>

[Bruno Chareyre (bruno-chareyre)]

>Christian should use a fixed time

I'd say GSTS should work even in that case (provided it is the source of the crash), and it is not too difficult to fix.
I'm not really sure PWave timestep is accurate for Hertz contact, is it?

[Chiara Modenese (chiara-modenese)]

On 28 November 2011 16:15, Chareyre <email address hidden> wrote:

> >Christian should use a fixed time
>
> I'd say GSTS should work even in that case (provided it is the source of
> the crash), and it is not too difficult to fix.
> I'm not really sure PWave timestep is accurate for Hertz contact, is it?
>
I have used it without any problem so far in QS condition (and I have also
run tests for different time steps to check that the solution is truly not
time step dependent.) Besides that, a very similar definition is used by
many other authors (e.g. Thornton) with non linear contact models (e.g.
Hertzian law).
Chiara

>
> --
> You received this bug notification because you are a member of Yade
> developers, which is the registrant for Yade.
> https://bugs.launchpad.net/bugs/897237
>
> Title:
> bug in neverErase?!
>
> Status in Yet Another Dynamic Engine:
> New
>
> Bug description:
> If I set neverErase=True Yade crashes with following output:
>
> O.run(3000,True)
> RuntimeError: std::exception
>
> and also:
>
> Traceback (most recent call last):
> File "/home/me/YADE/YADE2970/lib/yade-bzr2970/py/yade/qt/__init__.py",
> line 192, in refreshEvent
> 15:01:13
> self.refreshValues()
>
> 15:01:30
> File "/home/me/YADE/YADE2970/lib/yade-bzr2970/py/yade/qt/__init__.py",
> line 265, in refreshValues
> 15:01:34
> s=int(t); ms=int(t*1000)%1000; us=int(t*1000000)%1000;
> ns=int(t*1000000000)%1000
> ValueError: cannot convert float NaN to integer
>
>
> You can reproduce the bug with the attached script (similar to the one i
> submitted with version 2960).
> First try to set Law2_ScGeom_MindlinPhys_Mindlin(neverErase=False), this
> will work.
> Then set neverErase=True and yade crashes.
>
> To manage notifications about this bug go to:
> https://bugs.launchpad.net/yade/+bug/897237/+subscriptions
>
> _______________________________________________
> Mailing list: https://launchpad.net/~yade-dev
> Post to : <email address hidden>
> Unsubscribe : https://launchpad.net/~yade-dev
> More help : https://help.launchpad.net/ListHelp
>

[Christian Jakob (jakob-ifgt)]

While using O.dt=0.5*utils.PWaveTimeStep() instead of GSTS I get

1681 FATAL yade.ThreadRunner /home/me/YADE/yade/core/ThreadRunner.cpp:31 run: Exception occured:
Body #0 has velocity==NaN!
1719 ERROR yade.qt4 py/wrapper/yadeWrapper.cpp:368 wait: Simulation error encountered.
Traceback (most recent call last):
  File "/home/me/YADE/YADE2970/bin/yade-bzr2970", line 183, in runScript
    execfile(script,globals())
  File "macropore_neverErase.py", line 114, in <module>
    O.run(1000,True)
RuntimeError: std::exception

I think it should work for both time step estimations, so I will try to solve it like Bruno said: "Hertz Ip2 should at least assign kn=ks=0 in that case."

[Bruno Chareyre (bruno-chareyre)]

Chiara, it is not too difficult to define a time step resulting in a stable integration, and as soon as it is stable it is usually time-step independant. What I meant by "accurate" was: as close as possible to the stability condition. What pre-factor "p" are you using in dt=p*PWave?

[Chiara Modenese (chiara-modenese)]

On 29 November 2011 08:38, Chareyre <email address hidden> wrote:

> Chiara, it is not too difficult to define a time step resulting in a
> stable integration, and as soon as it is stable it is usually time-step
> independant.

I agree, it was you to ask the question :-)

> What I meant by "accurate" was: as close as possible to the
> stability condition.

What do you mean by stability? What definition you are referring to now?
Thanks, Chiara

> What pre-factor "p" are you using in dt=p*PWave?
>
> --
> You received this bug notification because you are a member of Yade
> developers, which is the registrant for Yade.
> https://bugs.launchpad.net/bugs/897237
>
> Title:
> bug in neverErase?!
>
> Status in Yet Another Dynamic Engine:
> New
>
> Bug description:
> If I set neverErase=True Yade crashes with following output:
>
> O.run(3000,True)
> RuntimeError: std::exception
>
> and also:
>
> Traceback (most recent call last):
> File "/home/me/YADE/YADE2970/lib/yade-bzr2970/py/yade/qt/__init__.py",
> line 192, in refreshEvent
> 15:01:13
> self.refreshValues()
>
> 15:01:30
> File "/home/me/YADE/YADE2970/lib/yade-bzr2970/py/yade/qt/__init__.py",
> line 265, in refreshValues
> 15:01:34
> s=int(t); ms=int(t*1000)%1000; us=int(t*1000000)%1000;
> ns=int(t*1000000000)%1000
> ValueError: cannot convert float NaN to integer
>
>
> You can reproduce the bug with the attached script (similar to the one i
> submitted with version 2960).
> First try to set Law2_ScGeom_MindlinPhys_Mindlin(neverErase=False), this
> will work.
> Then set neverErase=True and yade crashes.
>
> To manage notifications about this bug go to:
> https://bugs.launchpad.net/yade/+bug/897237/+subscriptions
>
> _______________________________________________
> Mailing list: https://launchpad.net/~yade-dev
> Post to : <email address hidden>
> Unsubscribe : https://launchpad.net/~yade-dev
> More help : https://help.launchpad.net/ListHelp
>

[Christian Jakob (jakob-ifgt)]

I tried to fix the kn and ks:

 if (uN<0) {
  if (neverErase) {
    phys->kn = 0.0;
    phys->ks = 0.0;
    phys->shearForce = phys->normalForce = Vector3r::Zero();
  }
  else {scene->interactions->requestErase(id1,id2); return;}
 }

Now my script is running to the end, but it is still crashing after the last step.

ValueError: cannot convert float NaN to integer

I think there are more values, that should be set to zero.
What about kno, kso, kr, ktw and so on?

In the case with Pwavetimestep it is still crashing at the first step:

1682 FATAL yade.ThreadRunner /home/me/YADE/yade/core/ThreadRunner.cpp:31 run: Exception occured:
Body #0 has velocity==NaN!
1719 ERROR yade.qt4 py/wrapper/yadeWrapper.cpp:368 wait: Simulation error encountered.
Traceback (most recent call last):
  File "/home/me/YADE/YADE-testing/bin/yade-bzr2970", line 183, in runScript
    execfile(script,globals())
  File "macropore_neverErase.py", line 114, in <module>
    O.run(1000,True)
RuntimeError: std::exception

[Chiara Modenese (chiara-modenese)]

On 29 November 2011 08:52, Chiara Modenese <email address hidden> wrote:

>
>
> On 29 November 2011 08:38, Chareyre <email address hidden> wrote:
>
>> Chiara, it is not too difficult to define a time step resulting in a
>> stable integration, and as soon as it is stable it is usually time-step
>> independant.
>
> I agree, it was you to ask the question :-)
>
>
>> What I meant by "accurate" was: as close as possible to the
>> stability condition.
>
> Say that you decrease the time step of a benchmark case (defined using the
PWave time step with a factor of 0.3) by 1 or 2 or 3 orders of magnitude.
It the mechanical response of the system does not change then I would say
that the solution is stable. Is there any other verification you would
suggest to do to verify the accuracy of the calculation?

> What do you mean by stability? What definition you are referring to now?
> Thanks, Chiara
>
>
>
>> What pre-factor "p" are you using in dt=p*PWave?
>>
>> --
>> You received this bug notification because you are a member of Yade
>> developers, which is the registrant for Yade.
>> https://bugs.launchpad.net/bugs/897237
>>
>> Title:
>> bug in neverErase?!
>>
>> Status in Yet Another Dynamic Engine:
>> New
>>
>> Bug description:
>> If I set neverErase=True Yade crashes with following output:
>>
>> O.run(3000,True)
>> RuntimeError: std::exception
>>
>> and also:
>>
>> Traceback (most recent call last):
>> File "/home/me/YADE/YADE2970/lib/yade-bzr2970/py/yade/qt/__init__.py",
>> line 192, in refreshEvent
>> 15:01:13
>> self.refreshValues()
>>
>> 15:01:30
>> File "/home/me/YADE/YADE2970/lib/yade-bzr2970/py/yade/qt/__init__.py",
>> line 265, in refreshValues
>> 15:01:34
>> s=int(t); ms=int(t*1000)%1000; us=int(t*1000000)%1000;
>> ns=int(t*1000000000)%1000
>> ValueError: cannot convert float NaN to integer
>>
>>
>> You can reproduce the bug with the attached script (similar to the one i
>> submitted with version 2960).
>> First try to set Law2_ScGeom_MindlinPhys_Mindlin(neverErase=False), this
>> will work.
>> Then set neverErase=True and yade crashes.
>>
>> To manage notifications about this bug go to:
>> https://bugs.launchpad.net/yade/+bug/897237/+subscriptions
>>
>> _______________________________________________
>> Mailing list: https://launchpad.net/~yade-dev
>> Post to : <email address hidden>
>> Unsubscribe : https://launchpad.net/~yade-dev
>> More help : https://help.launchpad.net/ListHelp
>>
>
>

[Bruno Chareyre (bruno-chareyre)]

Christian,
Yes. GSTS is fine after with kn/ks, but there are also junk values comming from the contact law itself apparently (after #7), maybe because of the other stiffness values.

Chiara,
Stable = !unstable ;)
If you don't see an explosion and NaN values in a few timesteps, then it is stable. There is a very narrow range of dt value where it is not exploding while not completely stable, something like dt +/- 10%.
So, if you get NaN velocities for dt>=T, set dt=0.8*T and you are on the safe side.
If you have, let's say, dt=0.01*T, it is stable too but you are wasting cpu time.

[Anton Gladky (gladky-anton)]

Christian, you script does not crash, if you keep the radii's constant.
I think the more proper way is to erase the particle and create a new one
in the same position but with the different radius.

I have slightly modified your script, which does it.
Please, see an attachment.
It does not crash, but gives a couple of warnings. Not sure it is 100% correct,
but it works.

Please, have a look at it.

PS Sorry, I am not able to participate actively in discussion now.

Anton

[Christian Jakob (jakob-ifgt)]

thank you very much anton, i get some error messages, but it works.

3945 ERROR yade.InteractionContainer /home/me/YADE/yade/core/InteractionContainer.cpp:55 erase: InteractionContainer::erase: attempt to delete non-existent interaction ##43+380
3945 ERROR yade.InteractionContainer /home/me/YADE/yade/core/InteractionContainer.cpp:55 erase: InteractionContainer::erase: attempt to delete non-existent interaction ##126+134
3945 ERROR yade.InteractionContainer /home/me/YADE/yade/core/InteractionContainer.cpp:55 erase: InteractionContainer::erase: attempt to delete non-existent interaction ##131+300
3945 ERROR yade.InteractionContainer /home/me/YADE/yade/core/InteractionContainer.cpp:55 erase: InteractionContainer::erase: attempt to delete non-existent interaction ##316+367
3945 ERROR yade.InteractionContainer /home/me/YADE/yade/core/InteractionContainer.cpp:55 erase: InteractionContainer::erase: attempt to delete non-existent interaction ##326+372

P.S. it does not work with O.dt=0.5*utils.PWaveTimeStep() instead of GSTS

[Chiara Modenese (chiara-modenese)]

On 29 November 2011 09:14, Chareyre <email address hidden> wrote:

> Christian,
> Yes. GSTS is fine after with kn/ks, but there are also junk values comming
> from the contact law itself apparently (after #7), maybe because of the
> other stiffness values.
>
> Chiara,
> Stable = !unstable ;)
> If you don't see an explosion and NaN values in a few timesteps, then it
> is stable.

Then according to my tests is stable while using PWave time step. But then
I would not talk about "accuracy" only because there is no explosion. Hope
this answers your question.

> There is a very narrow range of dt value where it is not exploding while
> not completely stable, something like dt +/- 10%.
> So, if you get NaN velocities for dt>=T, set dt=0.8*T and you are on the
> safe side.
> If you have, let's say, dt=0.01*T, it is stable too but you are wasting
> cpu time.
>
> --
> You received this bug notification because you are a member of Yade
> developers, which is the registrant for Yade.
> https://bugs.launchpad.net/bugs/897237
>
> Title:
> bug in neverErase?!
>
> Status in Yet Another Dynamic Engine:
> New
>
> Bug description:
> If I set neverErase=True Yade crashes with following output:
>
> O.run(3000,True)
> RuntimeError: std::exception
>
> and also:
>
> Traceback (most recent call last):
> File "/home/me/YADE/YADE2970/lib/yade-bzr2970/py/yade/qt/__init__.py",
> line 192, in refreshEvent
> 15:01:13
> self.refreshValues()
>
> 15:01:30
> File "/home/me/YADE/YADE2970/lib/yade-bzr2970/py/yade/qt/__init__.py",
> line 265, in refreshValues
> 15:01:34
> s=int(t); ms=int(t*1000)%1000; us=int(t*1000000)%1000;
> ns=int(t*1000000000)%1000
> ValueError: cannot convert float NaN to integer
>
>
> You can reproduce the bug with the attached script (similar to the one i
> submitted with version 2960).
> First try to set Law2_ScGeom_MindlinPhys_Mindlin(neverErase=False), this
> will work.
> Then set neverErase=True and yade crashes.
>
> To manage notifications about this bug go to:
> https://bugs.launchpad.net/yade/+bug/897237/+subscriptions
>
> _______________________________________________
> Mailing list: https://launchpad.net/~yade-dev
> Post to : <email address hidden>
> Unsubscribe : https://launchpad.net/~yade-dev
> More help : https://help.launchpad.net/ListHelp
>

[Bruno Chareyre (bruno-chareyre)]

"Using PWave time step" doesn't tell what dt you are using.
Do you have dt=PWave, dt=0.1*PWave,... ?
In most example scripts, there is something like 0.5*PWave.

An "accurate" estimate of the critical timestep is something like T'=0.99*T, with T the theoretical stability value (defined here https://www.yade-dem.org/doc/formulation.html#dem-simulations). I was referring to accuracy in this sense.
But actually, as said above, getting a stable integration is also sufficient condition for getting an accurate integration in QS cases.

[Jérôme Duriez (jduriez)]

"T the theoretical stability value "

I was wondering if such values exists (I mean, can be theoretically derived) for any cases (more than two particles along one axis especially) ? The (your ?) paragraph you quoted Bruno made me consider the question...

I just realized that the idea is to consider the discrete numerical model with a vector X = ( x(body1),y(body1),z(body1),x(body2),y(body2),z(body2),.... ), and write that our DEM computations leads to write
             X(t+dt) = X(t) * f(dt) (1)
     With f(t) determined by the fact that the second derivative of X, d²X/dt² is given by a linear system M d²X/dt² + K X = 0, with M and K the "generalized mass and stiffness matrices".
    Then you end up saying that equation (1) is stable for dt values given by the classical relation.

Did I understand well ? If yes, I have still few questions :
* does this not depend on the integration scheme which is used ? Depending on how we compute X(t+dt) from X(t) and d²X/dt², we do not end up with same f(dt) in equation (1), isn't it ? Then, different values of critical values of dt ? Or not ?

* the consideration of a linear system M d²X/dt² + KX = 0 assumes no boundary conditions and no gravity, isn't it ?

* in any case, there should be a typo at last line of https://www.yade-dem.org/doc/formulation.html#general-mass-spring-system : m / k does not have the same unit as eigenvalue of M-1 * K
(it would be difficult for me to perform myself the correction)

[Christian Jakob (jakob-ifgt)]

I tried another script without facets (using aabbWalls) and without changing particle radii,
it crashes after 5 or 6 steps (when neverErase=True):

5468 FATAL yade.ThreadRunner /home/me/YADE/yade/core/ThreadRunner.cpp:31 run: Exception occured:
Body #3084 has velocity==NaN!

I think the problem is the changes, that Bruno implemented with neverErase. All contacts remain active with this option. But we need a rupture distance, at which liquid bridges break and the contact gets inactive. Something like

if (uN<0) && (uN>rupture_distance){
 if (neverErase) {
   phys->kn = 0.0;
   phys->ks = 0.0;
   #...maybe more content here...
   phys->shearForce = phys->normalForce = Vector3r::Zero();
 }
 else {scene->interactions->requestErase(id1,id2); return;}
}
if (uN<rupture_distance) {scene->interactions->requestErase(id1,id2); return;}

[Christian Jakob (jakob-ifgt)]

... This leads to another question:

Where is the calculation of the rupture distance? I could not find it in Law2_ScGeom_CapillaryPhys_Capillarity.cpp.
Does this mean, that the capillary force is calculated only for real contacts (when uN>0)?
It should be created even before particles come in contact (creation distance) and it should be erased, when particles achieve a critical distance (quasi-static rupture distance).

[Bruno Chareyre (bruno-chareyre)]

@Jérôme
Yes, I typed this paragraph (some parts inherit some Vaclav's text to).
1) Of course, critical is relative to a specific integration scheme. As you know implicit schemes don't even have a critical timestep (they have more accuracy problems however).
2) BCs don't change anything in the equations (if boundaries are dynamic they are ordinary bodies, if they are not dynamic they are out of the scope of this paragraph). Gravity does not change eigenfrequencies.
3) It is M.K-1 indeed, thanks!

@Christian

Rupture is (in our model) when there is no more mathematical solution to Laplace equation. There is no constant rupture distance, rupture is a function of D1/D2, distance, and succion. Practicaly, if the law result in Fcap=0, it means we are out of the range of acceptable parameters, then the bridge is broken.

For the rest, you are right. I have to check the code to answer more precisely. Be right back, I can't browse right now.

[Bruno Chareyre (bruno-chareyre)]

@Christian

The answer to #18 is at Law2_ScGeom_CapillaryPhys_Capillarity.cpp:163, there is no need to erase interactions in the contact law because the capillary law takes care of that:

if (!Vinterpol){
   ...
   scene->interactions->requestErase(id1,id2);
}

I.e. if the volume of liquid is zero, Laplace problem has no solution -> the bridge is broken -> we delete the interaction.
All good. This has been well tested in the last years.

The problem is we never tried to delete bodies, and at line 90 you can find the recent addition (I appended the FIXME):
if (b1 and b2) {// FIXME: possible memory leak here?

This test was initialy supposed to fix the crash after erasing bodies, but now if !b1 or !b2, the existing distant interactions will never be erased.

[Bruno Chareyre (bruno-chareyre)]

Oh wait... erasing bodies triggers erasing interactions if I'm not wrong, so there is still another problem.

[Bruno Chareyre (bruno-chareyre)]

Heh no, erasing bodies "requests" erasing interactions, but it is delayed by collider stride, so I think it is the problem. I'll try and fix.
Chrsitian, I downloaded your script above, I don't understand why hertzOn=0 while hertz is used.

[Anton Gladky (gladky-anton)]

> Heh no, erasing bodies "requests" erasing interactions, but it is delayed by collider stride, so I think it is the problem. I'll try and fix.

I do not see any problems there.
I think, the warning appears, when both interacted bodies are deleted
at the same step.

IMHO, no reasons to break a working code.

Anton

[Bruno Chareyre (bruno-chareyre)]

>no reasons to break a working code.

I don't think I will touch the part of the code you are afraid I could touch. ;)

[Bruno Chareyre (bruno-chareyre)]

Christian, I tried to delete bodies with capillary interactions using a fresh script and I see no problems after bzr2973.
Check what's inside the attached script, there a few things that should also be included in yours as explained in the documentation:

"In order to allow capillary forces between distant spheres, it is necessary to enlarge the bounding boxes using :yref:`Bo1_Sphere_Aabb::aabbEnlargeFactor` and make the Ig2 define define distant interactions via :yref:`interactionDetectionFactor<Ig2_Sphere_Sphere_ScGeom::interactionDetectionFactor>`. It is also necessary to disable interactions removal by the constitutive law (:yref:`Law2<Law2_ScGeom_FrictPhys_CundallStrack::neverErase>=True`)."

I also included automatic definition of hertzOn, so it is not necessary to define it (if you define it, it will be overwritten anyway).
You can try the script attached above and turn neverErase=1, then you will get distant interactions (see also createDistantMeniscii flag)

[Christian Jakob (jakob-ifgt)]

Thanks Bruno, my script is now working without crash.