# I calculated the kinetic energy by E=0.5mv^2+0.5Iw^2 but it is not equal to kineticEnergy() when O.cell.velGrad is set

I am conducting a simulation where the energy results are important for the study. So I want to study how the kinetic energy is calculated in YADE. In the documentation it shows kinetic energy is calculated by E=0.5*m*v^2+0.5*I*w^2, which is definitely normal and accurate. So I have tried to calculated kinetic energy by myself, as the particle mass m, velocity v, inertia I and angular velocity w can be all obtain by using O.bodies[]state.mass, O.bodies[].state.vel, O.bodies[].state.inertia and O.bodies[].state.angVel. But when I have set the value of O.cell.velGrad, the E=0.5*m*v^2+0.5*I*w^2 calculated by myself is much larger than utils.kineticEnergy(). And when I didn't use O.cell.velGrad, the E=0.5*m*v^2+0.5*I*w^2 calculated by myself is the same with utils.kineticEnergy()

I am thinking why this happened? Why setting the value of O.cell.velGrad will make E=0.5*m*v^2+0.5*I*w^2 larger than the utils.kineticEnergy()?

My code for calculating the E=0.5*m*v^2+0.5*I*w^2 is:

kineE=0
for i in bodies_id:
m=O.bodies[i].state.mass
v2=numpy.dot(O.bodies[i].state.vel, O.bodies[i].state.vel)
I=O.bodies[i].state.inertia
w2=numpy.dot(O.bodies[i].state.angVel, O.bodies[i].state.angVel)
kineE=kineE+0.5*m*v2+0.5*I*w2

Thank you very much!!
Yuxuan

## Question information

Language:
English Edit question
Status:
Solved
For:
Assignee:
No assignee Edit question
Solved by:
Yuxuan Wen
Solved:
Last query:
 Revision history for this message Robert Caulk (rcaulk) said on 2021-08-16: #1

Hello,

Cheers,

Robert

 Revision history for this message Yuxuan Wen (wenyuxuan) said on 2021-08-16: #2

Yes, I set the velGrad in the format as shown in the documentation: O.cell.velGrad=Matrix3(-0.5,0,0,0,-0.5,0,0,0,-0.5). And then the kinetic energy calculated by E=0.5*m*v^2+0.5*I*w^2 is different with the kineticEnergy(). I made a simplified code of the simulation, as shown follows. When you run this code directly, the kinetic energy is different, and when you comment the line of O.cell.velGrad, the kinetic energy will be the same.

#################################
from yade import pack, plot, qt, export, os
O.periodic=True
O.cell.hSize=Matrix3(0.03,0,0, 0,0.03,0, 0,0,0.03)
sp=pack.SpherePack()
sp.makeCloud(Vector3(0,0,0), Vector3(0.03,0.03,0.03), num=1000, rMean=0.001, rRelFuzz=0, distributeMass=False, periodic=True, seed=1)

O.engines=[
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom(), Ig2_Box_Sphere_ScGeom()],
[Ip2_FrictMat_FrictMat_MindlinPhys(betan=0.3,betas=0.3)],
[Law2_ScGeom_MindlinPhys_Mindlin()] ),
NewtonIntegrator(gravity=(0,0,-9.81),damping=0.0),
PyRunner(command='check()',iterPeriod=1000),
]
O.dt=0.3*utils.PWaveTimeStep()
O.cell.velGrad=Matrix3(-0.5,0,0, 0,-0.5,0, 0,0,-0.5) # if comment this line, kineE will be the same with kineticEnergy()

def check():
kineE=0
for i in O.bodies:
m=i.state.mass
v2=numpy.dot(i.state.vel, i.state.vel)
I=sqrt(numpy.dot(i.state.inertia,i.state.inertia))
w2=numpy.dot(i.state.angVel, i.state.angVel)
kineE=kineE+0.5*m*v2+0.5*I*w2 # compare with utils.kineticEnergy()
print(kineE, kineticEnergy())
#################################

Thank you so much!

Kind Regards,
Yuxuan

 Revision history for this message Jan Stránský (honzik) said on 2021-08-17 #3

Hello,

both values are computed differently.
E=0.5*m*v^2+0.5*I*w^2
of course still holds, but the values of "v" and "w" are different.
Have a look at source code : conditions "if (scene->isPeriodic)", variables like "spin = scene->cell->getSpin()" and their usage, comments "Only take in account the fluctuation velocity" etc.
If still anything is not clear, please let us know.

If you think this discrepancy should be mentioned in the documentation, feel free to update it yourself or create an issue on gitlab pages of the project (not to be forgotten).

Cheers
Jan

 Revision history for this message Yuxuan Wen (wenyuxuan) said on 2021-08-17: #4

Hello Jan,

Thank you for you help and sending me the position of the source code! I am much appreciated!

But I am still a little confused with the definition of the fluctuation velocity. As you said, the E=0.5*m*v^2+0.5*I*w^2 still holds, and "v" and "w" are the fluctuation velocity . In the source code, the fluctuation velocity is defined as:
bodyFluctuationVel(state->pos - state->vel * scene->dt, state->vel, scene->cell->velGrad)
which is a function of 1)updated position, 2) v, and 3) velGrad.
But I don't know how to calculate the fluctuation velocity by these 3 terms. I searched the documentation by "bodyFluctuationVel" and didn't find anything matched. I wonder if it is possible that you could tell me the expression of the fluctuation velocity? Thank you!

Kind Regards,
Yuxuan

 Revision history for this message Yuxuan Wen (wenyuxuan) said on 2021-08-17: #5

Hello Jan,

Thanks for your suggestion, now I have created an issue on gitlab pages of the project . Because I am still confused with the definition of fluctuation velocity, so I didn't update it myself but created an issue on the gitlab pages.

Kind Regards,
Yuxuan

 Revision history for this message Jan Stránský (honzik) said on 2021-08-19: #6

Probably again only source code is relevant here .

The idea is this:
Let's have a particle at (0,0,0), i.e. one corner of periodic cell.
Properties of the particle should (must) be identical with its periodic images (in this case e.g. particles in other corners of the periodic cell). See  and pictures there for a reference.
Let's say the (0,0,0) particle has zero velocity.
Now consider the cell has some velGrad. E.g. it is stretching in x direction. Then he "image particles" in the "right face" of the cell has some non-zero positive velocity.
Now, from the periodicity point of view, it should not matter if you choose the (0,0,0) particle or any of its periodic images.
The (0,0,0) particle has zero velocity, the periodic images have non-zero velocity, which one should be chosen?? The "fluctuation part" of the velocity is chosen in Yade.

Similarly it is with angular velocity and cell spin 

Cheers
Jan

 Revision history for this message Yuxuan Wen (wenyuxuan) said on 2021-08-19: #7

Thank you so much Jan!! Now I understand the definition of fluctuation velocity and why it is chosen in Yade!
Thanks again for you help!! I am much appreciated!
Yuxuan