plastic energy dissipation

sorry for the mistakes. Here is the new version of my code:

aa=dict(O.energy.items())
E=aa['plastDissip']
f=open('plastic energy dissipation.txt','w')

for k in range (1,4000):
        O.run(1,True)
        aa=dict(O.energy.items())
        for i in O.interactions:
                Ft=i.phys.shearForce
                du=i.geom.shearInc
                kt=i.phys.ks
                Fn=i.phys.normalForce
                Ftt=Ft-kt*du
                maxFs=Fn.norm()*i.phys.tangensOfFrictionAngle

                if(Ftt.norm() > maxFs):
                        ratio = maxFs / Ftt.norm()
                        trialForce=Ftt
                        Ftt *= ratio
                        e=((1/kt)*(trialForce-Ftt)).dot(Ftt)
                        E=E+e

        f.write('%d %f %f %f\n' %(O.iter,O.time,E,aa['plastDissip']))

f.close()

Hello Lingran,

actually it is very lucky that they are different :-) In c++ code, the
values are transfered from previous time step to the current one, while in
your python script, your values are always the current one. For example, Ft
= i.phys.shearForce is NOT the same as the "equivalent" in c++ ( Vector3r&
shearForce = geom->rotate(phys->shearForce); ). even neglecting the
rotation. shearForce in c++ is one time step earlier than in the python.
So, for one interaction, the correct approach should be:

Ft=i.phys.shearForce
O.step() # now you use the same value of Ft as in c++
du=i.geom.shearInc
kt=i.phys.ks
Ftt=Ft-kt*du
maxFs=Fn.norm()*i.phys.tangensOfFrictionAngle
if(Ftt.norm() > maxFs):
  ratio = maxFs / Ftt.norm()
  trialForce=Ftt
  Ftt *= ratio
  e=((1/kt)*(trialForce-Ftt)).dot(Ftt)
print Ftt,i.phys.shearForce # now both should be the same

For more interactions, you can store shearForces of the interactions in a
list and after O.step() use it in for loop over interaction.

A few commets to your script (not very important):
- you can access O.energy['plastDissip'] directly (without need of
O.energy.items() )
- during the simulation, you can store data to plot.data (with plot.addData
function, e.g.
   plot.addData(i=O.iter,t=O.time,e=E,p=O.energy['plastDissip'])
and then save it with plot.saveDataTxt function. For non-alphabetical
ordered of variables, use O.saveDataTxt(vars=('i','t','e','p'))
- I would use O.step() insteat of O.run(1,True), but it really does not
matter :-)

cheers
Jan

2013/8/16 lingran <email address hidden>

> Question #234129 on Yade changed:
> https://answers.launchpad.net/yade/+question/234129
>
> lingran gave more information on the question:
> sorry for the mistakes. Here is the new version of my code:
>
> aa=dict(O.energy.items())
> E=aa['plastDissip']
> f=open('plastic energy dissipation.txt','w')
>
> for k in range (1,4000):
> O.run(1,True)
> aa=dict(O.energy.items())
> for i in O.interactions:
> Ft=i.phys.shearForce
> du=i.geom.shearInc
> kt=i.phys.ks
> Fn=i.phys.normalForce
> Ftt=Ft-kt*du
> maxFs=Fn.norm()*i.phys.tangensOfFrictionAngle
>
> if(Ftt.norm() > maxFs):
> ratio = maxFs / Ftt.norm()
> trialForce=Ftt
> Ftt *= ratio
> e=((1/kt)*(trialForce-Ftt)).dot(Ftt)
> E=E+e
>
> f.write('%d %f %f %f\n' %(O.iter,O.time,E,aa['plastDissip']))
>
> f.close()
>
> --
> 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
>

Thanks a lot, Jan, you are saving my life.
But I am still a little confused. I did a small test, it turns out i.phys.shearForce is equal to Ft rather than Ftt.
Did you miss something after O.step()?
see:

i=O.interactions[975,2702]
Ft=i.phys.shearForce # -- Ft=(-263.88, -780.43,0)
O.step() # now you use the same value of Ft as in c++
du=i.geom.shearInc
kt=i.phys.ks
Ftt=Ft-kt*du # --Ftt=(-274.83, -812.82, 0)
maxFs=Fn.norm()*i.phys.tangensOfFrictionAngle #--maxFs=223.55
If (Ftt.norm() > maxFs):
      ratio = maxFs / Ftt.norm()
      trialForce=Ftt
      Ftt *= ratio # Ftt=(-71.61, -211.77,0)
      e=((1/kt)*(trialForce-Ftt)).dot(Ftt)
      i.phys.shearForce # -- i.phys.shearForce =(-263.88, -780.43,0)

Best regards:
Lingran

Hmm, it seems I deleted i.phys.normalForce.. should be:

Ft=i.phys.shearForce
O.step()
du=i.geom.shearInc
kt=i.phys.ks
Fn=i.phys.normalForce
Ftt=Ft-kt*du
maxFs=Fn.norm()*i.phys.tangensOfFrictionAngle
if(Ftt.norm() > maxFs):
  ratio = maxFs / Ftt.norm()
  trialForce=Ftt
  Ftt *= ratio
  e=((1/kt)*(trialForce-Ftt)).dot(Ftt)
  E=E+e

> i=O.interactions[975,2702]
> Ft=i.phys.shearForce # -- Ft=(-263.88, -780.43,0)
> O.step() # now you use the same value of Ft as in c++
> du=i.geom.shearInc
> kt=i.phys.ks
> Ftt=Ft-kt*du # --Ftt=(-274.83, -812.82, 0)
> maxFs=Fn.norm()*i.phys.tangensOfFrictionAngle #--maxFs=223.55
> If (Ftt.norm() > maxFs):
> ratio = maxFs / Ftt.norm()
> trialForce=Ftt
> Ftt *= ratio # Ftt=(-71.61, -211.77,0)
> e=((1/kt)*(trialForce-Ftt)).dot(Ftt)
> i.phys.shearForce # -- i.phys.shearForce =(-263.88, -780.43,0)
>
>
the missing of Fn may be a problem aslo here. i.phys.shearForce.norm() >
maxFs, which should not be..

In case it does not help, please insert to the email the working script (if
it is too large, only the important parts, but such that it works) so I (or
others) can try it.
Thanks
Jan

Hi, Jan
Sorry for replying to you so late, it’s not possible to run yade simulations on my personal windows 8 computer at home. I should replace windows 8 with windows 7 as soon as possible.

Here is the whole script: http://pastebin.com/xM5Djdmq

Now I am still confused about how plastic energy dissipation is calculated in c++ code and python respectively?

1. “In c++ code, the values are transfered from previous time step to the current one, while in your python script, your values are always the current one”.

The question is I also use du=i.geom.shearInc and calculate Ft=Ft-kt*du, why it can not be considered as a updating of value?

2. I got same values of elastic energy at each O. iter both from c++ and python. If c++ is one time step earlier than python, then the values shouldn’t be the same.

3. It is still difficult for me to understand how does your script work?

        print Ftt,i.phys.shearForce # now both should be the same----- do you mean Ftt=i.phys.shearForce here? Then why should them be equal to each other?

Thanks.
Best regards.
Lingran

2013/8/19 lingran <email address hidden>:
> Sorry for replying to you so late, it’s not possible to run yade simulations on my personal windows 8 computer at home. I should replace windows 8 with windows 7 as soon as possible.

Yade runs on Windows???

Anton

Hi Lingran,

> Sorry for replying to you so late,

not problem at all for me :-)

> it’s not possible to run yade simulations on my personal windows 8
> computer at home. I should replace windows 8 with windows 7 as soon as
> possible.
>

you are able to run yade on windows 7 (just for my curiosity)?

>
> Here is the whole script: http://pastebin.com/xM5Djdmq

thanks, I will try it

>
>
> 1. “In c++ code, the values are transfered from previous time step
> to the current one, while in your python script, your values are always
> the current one”.
>
> The question is I also use du=i.geom.shearInc and calculate Ft=Ft-kt*du,
> why it can not be considered as a updating of value?
>

Basicly you are updating the value, you are right. But the starting point
of the update is different. In c++, you start with shearForce from previous
step and do some update (to current step). In python you start to update
already updated (in c++) shearForce (in current step). Maybe I got
something wrong, that's why I asked you for the script to test it myself,
but just now I see this difference as the most important one.

>
> 2. I got same values of elastic energy at each O. iter both from
> c++ and python. If c++ is one time step earlier than python, then the
> values shouldn’t be the same.
>

elastic energy is computed from current values of force (i.e. the same in
c++ and python), so this is ok

>
> 3. It is still difficult for me to understand how does your script
> work?
>
> print Ftt,i.phys.shearForce # now both should be the same-----
> do you mean Ftt=i.phys.shearForce here? Then why should them be equal
> to each other?
>
>
sorry, I will try your script and then try to formulate the answer in as
understandable way as possible :-)

cheers
Jan

Hello Lingran,
I tried your script, but it runs too long for playing.. Would it be
possible to shorten it (e.g. to use not so many particles and not so many
time steps), but such that is still reproduces your problem?
Thanks
Jan

Sorry for the long and complex script. Here is the short version: http://pastebin.com/TLVk040M.
Besides, for your curiosity, yade can’t be run on windows. I mentioned windows 7 just because it is easier to install Ubuntu on it compare to windows 8. Luckily, somebody is helping me replace windows 8 by windows 7, Ubuntu and yade are on the way to my laptop :-).

By proposing question here, I seek help to confirm energy balance, that means:
Ek + Es + Emgh + Edamping + Eplasitc=constant,
where Ek = kinetic energy, Es = elastic energy( including normal and shear components), Emgh = gravity potential energy, Edamping =non-viscous damping energy (can be ignored if we set damping=0), Eplastic = plastic energy dissipation.
I have no doubt of calculating Ek, Es, Emgh, Edamping, but the most difficult part is Eplastic. First, the value I calculated is not consistent with O.energy[‘PlastDissip’], second, I haven't find a good Eplastic value to balance energy, not yet…..

Thanks.
Best regards.
Lingran

Hi Lingran,

Even before considering plastic dissipation, do you get a correct energy balance when friction=0 (no dissipation)?
I guess you don't. The reason is twofold.

(1) the global energy tracking is present in Law2_ScGeom_FrictPhys_CundallStrack but not in Law2_ScGeom6D_CohFrictPhys_CohesionMoment (*). Since you are combining them, you see something in O.energy['elastPotential'] but it only partly reflects the total elastic energy.
You can get elastic energy independently for each contact type with lawN.normElastEnergy() and lawN.shearElastEnergy() (where N will be 1 and 2 in your script).

(2) you are using contact moments, and unfortunately there is no law2.bendingElastEnergy(). So this energy is simply missing in the total, even after considering (1) above. You may write this function in c++ actually, it would be useful for everyone. It is not very difficult using the example of normElastEnergy(). Let us know.

Bruno

(*) It was my decision to not track elastic energy since, unlike plastic dissipation, elastic energy does not have to be incremented at each step. Computing it all the time was significantly slowing down the simulations for no good.

Hi Bruno,
Thanks for your comments. You are right and I will try my best to add law2.bendingElastEnergy() in c++ code.

But I still have several doubts.
1. like Jan said, Why c++ is one time step earlier than python?
2. When use Law2_ScGeom6D_CohFrictPhys_CohesionMoment() and CohFrictMat, why do we need to set both alphaKr and alphaKtw to a non-zero value, otherwise the simulation would crush?
3. Energy is not balanced even without considering friction and moments, that means: utils.kineticEnergy+ Emgh(sum of (i.state.mass*(i.state.pos[1]*9.81 ))+ law2.normElastEnergy()+law2.shearElastEnergy()+O.energy[‘nonviscDamp’]!=constant.. The error increases with the increasing number of particles.

Thanks.
Best regards.
Lingran

1. I will only put Jan's comment a bit differently. The "du" that should be used in "Ftt=Ft-kt*du" is the future displacement between t and t+dt. Since t+dt is not computed yet, du is not know. So you are using du of the previous step, but it is inconsistent.

2. I didn't know but searching alphaKr in the code shows this:
Real AlphaKr = 2.0*sdec1->alphaKr*sdec2->alphaKr/(sdec1->alphaKr+sdec2->alphaKr)

division by zero = no good... A small bug to fix.
It is not really a problem for you. If you do nothing (keep all defaults) there will be no moment (momentRotationLaw=False) and it will not crash (default alpha's are not zero).

3. There are a few known sources of energy changes. For instance creating a contact leads to unbalanced energy change. How big is the error compared to total energy? Does it decrease with the timestep?
I know energy is conserved very accurately in dense quasi-static systems. Is it what you are interested in, or do you need also energy conservation in dynamic problems?

Hi Bruno,
Thanks for your answers. They do help me a lot.

I am interest in dynamic energy conservation issues. I am dealing with impact of a boulder onto a soil layer. In order to investigate impact mechanisms, the soil layer is divided into different parts. Energy transmission and dissipation inside each soil layer supposes to be tracked, and before doing that, the total energy should be balanced.

For your comments, I did several tests and found that the errors do decrease with time step. When simply conducting gravity deposit, a sample of 2000 particles, without friction and moments, the error is about 0.03%. But when dealing with impact problems, consider friction (use O.energy[‘plastDissip’]) and moments, the error sometimes could be as high as 3~10%, in addition, the error depends on the size of the whole sample, the number of particles, impact energy, time step, etc.

Thanks.
Lingran

Hi all,

As promised, I try to add bendingElastEnergy and twistElastEnergy to Law2_ScGeom6D_CohFrictPhys_CohesionMoment. As a result, the corresponding cpp and hpp files are slightly modified. In addition, the new source codes are recompiled to update yade.

Several simple 2D tests are conducted to check this extended function. The results show that when Law2_ScGeom6D_CohFrictPhys_CohesionMoment is only used, bendingElastEnergy and twistElastEnergy values are verified and consistent with calculated values, while in the other hand, when Law2_ScGeom6D_CohFrictPhys_CohesionMoment is combined with Law2_ScGeom_FrictPhys_CundallStrack, bendingElastEnergy is given as ‘nan’ or ‘inf’ in yade.

Here are the diff files:
cpp_diff_for_cohesiveLaw: http://pastebin.com/eWzpLdGC
hpp_diff_for_cohesiveLaw: http://pastebin.com/XD18Wktn

Here are the modified source codes:
cpp_add bending_to_cohesiveLaw: http://pastebin.com/XQADy5Vt
hpp_add bending_to_cohesiveLaw: http://pastebin.com/GuUiL2g0

Hope it helps.
Best regards.
Lingran

Wow! Thanks a tone. We shall include that.

You are lucky that it does not simply crash when combined with
Law2_ScGeom_FrictPhys_CundallStrack.
I think you are actually trying to get phys->moment_bending for an
object that does not have the member moment_bending (a FrictPhys object).

You need to filter the interactions. An efficient (i.e. avoiding
dynamic_casts) method is this:
if (phys->getClassIndex() == CohFrictPhys::getClassIndexStatic()) {
//... do something with CohFrictPhys object }

Another interesting addition would be
Law2_ScGeom6D_CohFrictPhys_CohesionMoment::elasticEnergy(), since
currently we are forced to compute the total in the scripts...
It could be also added to Law2_ScGeom_FrictPhys_CundallStrack, in fact,
then overloaded in CohesionMoment.

Would you please fix that before inclusion in trunk?

Cheers.

Bruno

Hi Bruno,

yes, you are right, we should filter interactions when calculate bending elastic energy.
I will try my best.

Best regards.
Lingran

> You need to filter the interactions. An efficient (i.e. avoiding
> dynamic_casts) method is this:
> if (phys->getClassIndex() == CohFrictPhys::getClassIndexStatic())

I am almost sure that dynamic_cast is just as fast, or likely faster, than getClassIndex() - it is a virtual method, so the exact class of the object has to be looked up anyway. Try and measure.

@Vaclav
Interesting, thanks.
Now it is more clear for me why functorCache is so necessary.
Does it mean there is no better way in the present situation?

Would it make sense to add an "int index" to each Indexable and
initialize it once, like:
serializable.index = serializable.getClassIndex(); //? (just speculating)

Bruno

Hello Jan and other users,

The following is the new code I used to calculated plastic energy dissipation.
The result is quite close to yade's O.energy['plastDissip'] value, but huge computation time is needed especially when used to a big sample of 200,000 particles.

Could you please have a look? there might be some mistakes inside, and it would be great if someone could improve it.

E=O.energy['plastDissip']
e=0
energy=open('energy.txt','w')
for k in range (0,100):
 a=open('state_1.txt','w')
 b=open('state_2.txt','w')
 for i in O.interactions:
    Ft=i.phys.shearForce
    a.write('%d %d %f %f %f\n' %(i.id1,i.id2,Ft[0],Ft[1],Ft[2]))

 a.close()

 O.step()

 for i in O.interactions:
  Ftt=i.phys.shearForce
  du=i.geom.shearInc
  kt=i.phys.ks
  Fn=i.phys.normalForce
  maxFs=Fn.norm()*i.phys.tangensOfFrictionAngle
  b.write('%d %d %f %f %f %f %f %f %f %f\n' %(i.id1,i.id2,Ftt[0],Ftt[1],Ftt[2],maxFs,du[0],du[1],du[2],kt))

 b.close()

 x = numpy.loadtxt('state_1.txt')
 y = numpy.loadtxt('state_2.txt')

 for i in range (0,len(y)):
  Ftt=Vector3(y[i][2],y[i][3],y[i][4])
  maxFs=y[i][5]
  du=Vector3(y[i][6],y[i][7],y[i][8])
  kt=y[i][9]
  Ft=Vector3(0,0,0)
  for j in range (0,len(x)):
   if y[i][0]==x[j][0] and y[i][1]==x[j][1]:
       Ft=Vector3(x[j][2],x[j][3],x[j][4])
  Ftc=Ft-kt*du
  if Ftc.norm()> maxFs:
       #Ftt=Ftt/Ftt.norm()*maxFs
       due=-(Ftt-Ft)/kt
       #e+=Ftt.norm()*(du.norm()-due.norm())
       e+=Ftt.norm()*(du-due).norm()
 energy.write('%f %f %f\n'%(O.time,e,O.energy['plastDissip']-E))

energy.close()

Thanks and regards.
lingran

Sorry for the unclear code, this one is much better: http://pastebin.com/t8HUf9VS

Thanks.
Lingran

Hi Lingran,
May I ask a dumb question: what are you trying to achieve with such
code? Is it a verification of O.energy['plastDissip']?
On which aspects should/could it be improved?

Bruno

ps. Where are we for the rotational elastic energy? Would you like to
include a final version?

Hi Bruno and other users,

Yes, The code is supposed to be verified by obtaining the same value of O.energy['plastDissip'],then this method will be used futher to calculate local energy dissipations corresponds to only certains parts of interactions. Because in yade, we can track O.energy['plastDissip'] very easily, but it's only a global vaule, while in my case, I need to divide the soil sample into different layers and then calculate local plastic energy dissipation inside each layer.

so the final aim is to obtain local plastic energy dissipation values. In my opinion(maybe I am wrong), to achieve this aim, we need to filter some interactions and only consider the residual ones. This can be achieved by two ways: modify c++ codes or progamme python codes. I tried to filter interactions in c++ code but didn't succeed. then I tried to programme some python code, it works but it needs huge computation time. I hope this python code be improved so that it needs much less computation time.

Modify c++ codes and filter interactions has a splendid future, I needs some c++ knowledges and will keep trying.
It would greate if someone could offer me some helps.

Hope my problem doesn't try you guys crazy:)
Have a nice day!

Regards.
Lingran

Thank you for clarification. Getting the energies in subdomains is indeed not directly possible with O.energy.
This is not crazy, don't worry. :)
For python, I'm afraid it will always take a lot more time than c++, so I guess you will have to deal with the c++.

I am not sure how to implement this practically, it also depends on your precise needs.
Do you plan to filter the interactions only on the basis of their location? One way would be to simply include a an option to the O.energy system, like a custom predicate that would be checked inside scene->energy->add().
All the other solutions I can imagine are boring because they would imply to duplicate the same changes in all constitutive laws (or would work for only one of them).

The problem is: if you want the dissipated energy in a different region, then you have to run the same simulation again from the very begining (same problem in your python script). I see no simple solution to this.