Update interaction physics after changing material of spheres

Hi,

This question already appeared some times (e.g. https://answers.launchpad.net/yade/+question/197661), and the only way is to looping manually over the interactions to recompute the stifnnesses (adpat suggestion of answer #5 of quoted link).

Jerome

Well, we can then think about adding "something global" to recompute
such parameters.

Anton

2013/11/7 jduriez <email address hidden>:
> Question #238797 on Yade changed:
> https://answers.launchpad.net/yade/+question/238797
>
> Status: Open => Answered
>
> jduriez proposed the following answer:
> Hi,
>
> This question already appeared some times (e.g.
> https://answers.launchpad.net/yade/+question/197661), and the only way
> is to looping manually over the interactions to recompute the
> stifnnesses (adpat suggestion of answer #5 of quoted link).
>
> Jerome
>
> --
> You received this question notification because you are a member of
> yade-users, which is an answer contact for Yade.
>
> _______________________________________________
> Mailing list: https://launchpad.net/~yade-users
> Post to : <email address hidden>
> Unsubscribe : https://launchpad.net/~yade-users
> More help : https://help.launchpad.net/ListHelp

What do you think of changing the IP2-Functor?
instead of:
        if(interaction->phys) return; // no updates of an already existing contact necessary
        shared_ptr<MindlinPhys> contactPhysics(new MindlinPhys());
        interaction->phys = contactPhysics;
        FrictMat* mat1 = YADE_CAST<FrictMat*>(b1.get());
        FrictMat* mat2 = YADE_CAST<FrictMat*>(b2.get());

something like:
        FrictMat* mat1 = YADE_CAST<FrictMat*>(b1.get());
        FrictMat* mat2 = YADE_CAST<FrictMat*>(b2.get());
        if(interaction->phys->mat1==mat1 and interaction->phys->mat2=mat2) return;
        shared_ptr<MindlinPhys> contactPhysics(new MindlinPhys());
        interaction->phys = contactPhysics;
        // update reference to materials of bodies in contact
        interaction->phys->mat1=mat1
        interaction->phys->mat1=mat1

This way a reference to the materials from the last timestep of the bodies in interaction has to be stored in the contactPhysics.

Maybe

utils.calm() # not sure if necessary
O.interactions.clear()
O.step()

could work?

Jan

2013/11/7 Anton Gladky <email address hidden>

> Question #238797 on Yade changed:
> https://answers.launchpad.net/yade/+question/238797
>
> Anton Gladky proposed the following answer:
> Well, we can then think about adding "something global" to recompute
> such parameters.
>
> Anton
>
>
> 2013/11/7 jduriez <email address hidden>:
> > Question #238797 on Yade changed:
> > https://answers.launchpad.net/yade/+question/238797
> >
> > Status: Open => Answered
> >
> > jduriez proposed the following answer:
> > Hi,
> >
> > This question already appeared some times (e.g.
> > https://answers.launchpad.net/yade/+question/197661), and the only way
> > is to looping manually over the interactions to recompute the
> > stifnnesses (adpat suggestion of answer #5 of quoted link).
> >
> > Jerome
> >
> > --
> > You received this question notification because you are a member of
> > yade-users, which is an answer contact for Yade.
> >
> > _______________________________________________
> > Mailing list: https://launchpad.net/~yade-users
> > Post to : <email address hidden>
> > Unsubscribe : https://launchpad.net/~yade-users
> > More help : https://help.launchpad.net/ListHelp
>
> --
> You received this question notification because you are a member of
> yade-users, which is an answer contact for Yade.
>
> _______________________________________________
> Mailing list: https://launchpad.net/~yade-users
> Post to : <email address hidden>
> Unsubscribe : https://launchpad.net/~yade-users
> More help : https://help.launchpad.net/ListHelp
>

@Alexander

At the moment interaction->phys->mat1/2 does not exist, isnt't it ?

Good idea Jan.
It results in an update of interaction physics. But the problem is that my sphere packing explodes after rerunning simulation.
I think that this is what we can expect since the stiffnesses "explode", too. They get multiplyed by a factor of 1000. This means that contact forces explode and spheres are going crazy.
Maybe it is not possible to fill up my bowl with the soft particles...

@jduriez: yes you are right. But following my last comment, it would not be rational to implement it, because it causes explosion or at least for smaller differences in stiffnesses a lot of shaking (at least I think so)

2013/11/7 Alexander Eulitz [Eugen] <email address hidden>:
> Maybe it is not possible to fill up my bowl with the soft particles...

Soft particles are giving much more large overlaps, as their hard
analogues in equilibrium state. That is why
you have an explosion, when the stiffness is getting higher.

Anton

You can avoid explosions by using calm() function, which sets velocities of all bodies to zero.
In addition with a pyrunner, it can help you ... You can try something like this:

O.engines=O.engines+[PyRunner(iterPeriod=1,command='calm()',label='calmRunner')]
O.run(1000,True)

calmRunner.iterPeriod = 9
O.run(5000,True)

calmRunner.iterPeriod = 9999
O.run(10000,True)

calmRunner.iterPeriod = 99999
O.run(20000,True)
calmRunner.dead=True #deactivate calm runner

What about increasing the stiffness incrementally, e.g. increase it by a
factor of 2 or 10 (you probably have to try), run a couple of time steps
with high numerical damping and repeat this step until you reach your final
value.

Another way of avoiding the explosion is to downscale your particle radius
so that the particles are not in contact any more, i.e. ones you reach
equilibrium with your low stiffness you block all DoFs and then you scale
down the radius based on real interactions.

HTH
Klaus

On Thu, Nov 7, 2013 at 8:56 PM, Christian Jakob <
<email address hidden>> wrote:

> Question #238797 on Yade changed:
> https://answers.launchpad.net/yade/+question/238797
>
> Christian Jakob proposed the following answer:
> You can avoid explosions by using calm() function, which sets velocities
> of all bodies to zero.
> In addition with a pyrunner, it can help you ... You can try something
> like this:
>
>
> O.engines=O.engines+[PyRunner(iterPeriod=1,command='calm()',label='calmRunner')]
> O.run(1000,True)
>
> calmRunner.iterPeriod = 9
> O.run(5000,True)
>
> calmRunner.iterPeriod = 9999
> O.run(10000,True)
>
> calmRunner.iterPeriod = 99999
> O.run(20000,True)
> calmRunner.dead=True #deactivate calm runner
>
> --
> You received this question notification because you are a member of
> yade-users, which is an answer contact for Yade.
>
> _______________________________________________
> Mailing list: https://launchpad.net/~yade-users
> Post to : <email address hidden>
> Unsubscribe : https://launchpad.net/~yade-users
> More help : https://help.launchpad.net/ListHelp
>

I did it this way now:
There is a PyRunner which calls settleRealPacking() periodically.

def settleRealPacking():
    O.interactions.clear() #update contac physics after material was changed
    newton.damping = 0.9 #introduce high numerical damping
    print "unbF:", utils.unbalancedForce() #is just producing NaN?
    utils.calm() #delete relative velocities
    if utils.unbalancedForce()<0.1 : #threshold for settlement
        newton.damping = dmp #reset damping to original value
        checker.command = 'startProcessing()' #go on with simulated process

It looks ok, but I have to check whether the result is like the settled packing that would arise from a simulation with the real material from beginning on.

Thank you,
Alex

I support Klaus suggestion to increase stiffness more progressively, this is the best way to preserve the history. Deleting interactions make you loose shear forces, then the resulting packing is not equilibrated.

Since everything is proportional, you can update contact properties yourself:

mat1.Young *= stiffer
for i in O.interactions:
   i.phys.kn*= stiffer

We could do that in the Ip2 functor automatically as you suggest, but it doesn't need to add references (besides, your test would not work. interaction->phys->mat1 is still equal to mat1 even after changing a value in mat1).

We could overload the "write" operator on materials parameters so that a flag isChanged is turned on in the material object.
Then in Ip2 we would have something like:
if(interaction->phys && !mat1.isChanged && !mat2.isChanged) return;

A difficulty is to do that automatically for all variables of all material types. I don't see how we could do that. We should have in mind that resetting the interaction each time a variable is modified is not an option, in most case we want to change a parameter without deleting the existing contacts.

Last thing, in your case you can prevent explosions by increasing gravity with the same factor. It can be turned back to its initial value at the end (but again: progressively).

Bruno

Ok, thanks for the input. I see that this approach will be more correct then tablua rasa.
But will this progressive change of contact stiffness really be faster then using real material right from beginning?

I think that in case of HertzMindlin constitutive law, it should be
    i.phys.kno*= stiffer
Since contact normal force is calculated from kno and not kn, as it is with Cundall's constitutive law.
My test could maybe work if the 'equals' operator for FrictMats a and b would be overloaded with comparing all material parameters of a and b.

>will this progressive change of contact stiffness really be faster then using real material right from beginning?

Yes, definitely, if it is used the right way. Typically, you want to increase stiffness each 10 steps or so. If you increase at each step the loop for updating contact properties could take a significant time.

Yet, it is maybe not the best/simplest method. Actually, the computation time in your gravity deposition is ruled by gravity/stiffness. So, it would probably give the same benefit to increase gravity during deposition than to decrease stiffness.
Like for stiffness, you have to release gravity progressively in order to avoid explosions. The difference is that you don't have to take care of materials/interactions.
newton.gravity *= 0.99 #simple enough

ok, thanks Bruno.
I implemented it using the gracity approach and it works.