# save the evolution of the average kinetic energy of the grains and their average distance

hello, i wrote a code python of periodic inclined channel and i want to save the evolution of the average kinetic energy of the grains and their average distance every 500 period. when i open the file i found that kineticenergy=2 ( maximum value )during all simulation. this is my code:
#!/usr/bin/python
#Import libraries

#Inclination angle of the channel slope
slope=35

#Evolution of the kinetic energies of the flow over the average distance traveled by the grains.
n=7200 # number of spheres

# interactions
O.engines=O.engines+[NewtonIntegrator(damping=0.2, gravity=Vector3(-1*sin(slope),0.0,-1*cos(slope)), label='newtonIntegr'),
PyRunner(iterPeriod=500,command='AddData()'),#save every 500 periode of simulation
PyRunner(iterPeriod=23607,command='Data()'),# end of simulation
]

""" open file"""
fil=open('inclined_channel_40w_.dat', 'w')
fil.write('#Ec\t x \n')

x_g=0
for b in O.bodies:
x_g=x_g+b.state.pos
fil.write('%d %.3f\n' % (e,x_g/n))

#saving simulation
def Data():

#Timestep
O.dt=.1e-2
O.run(23707)
qt.View()
#utils.waitIfBatch()

## Question information

Language:
English Edit question
Status:
Solved
For:
Assignee:
No assignee Edit question
Solved by:
Jan Stránský
Solved:
2020-01-25
Last query:
2020-01-25
2020-01-20
 Jan Stránský (honzik) said on 2020-01-19: #1

Hello,

please provide the whole code you are using , without it we know nothing about your simulation and cannot help..

cheers
Jan

 hafsa (sebbah.hafsa) said on 2020-01-19: #2

hello Jan,
I divided my problem in two codes, the first is an initialization of configuration periodic channel without inclination:
#!/usr/bin/python
#Import libraries
from yade import pack, qt ,utils ,plot

# Material of Particles and walls

O.materials.append(ViscElMat(cn=33.5, cs=16.75, density=1.91, frictionAngle=0.4637, kn=2e5, ks=1e5, label='Mat'))

#Configuration:
length=40 #Streamwise length of the periodic cell, in diameter
width=40 #Spanwise length of the periodic cell, in diameter
height=2*length #heigth of the periodic cell
gravityVector = Vector3(0,0,-1) #Gravity vector

#Definition of the semi-periodic cell
O.periodic = True
O.cell.setBox(length,width,height)

# Reference walls: build two planes
wall_1=O.bodies.append(box(center=(20,0.1,9.5),extents=(20,0,9),fixed=True,wire=False,color=(1,0,1),material='Mat'))
wall_2=O.bodies.append(box(center=(20,39.9,9.5),extents=(20,0,9),fixed=True,wire=False,color=(1,0,1),material='Mat'))

# Build bumpy base made by particles
for i in range(1,81,1):
for j in range(2,82,1):
j=+1
i=+1

#Create a loose cloud of particle inside the cell
partCloud = pack.SpherePack()
partCloud.makeCloud(minCorner=(2,2,14),maxCorner=(length-5,width-5,height),rRelFuzz=0.,rMean=r_Prt,num =7200)
partCloud.toSimulation(material='Mat') #Send this packing to simulation with material Mat

# Simulation LOOP
O.engines = [
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Wall_Aabb(), Bo1_Facet_Aabb(),Bo1_Box_Aabb()],label='contactDetection',allowBiggerThanPeriod = True),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom(), Ig2_Box_Sphere_ScGeom()],
[Ip2_ViscElMat_ViscElMat_ViscElPhys()],
[Law2_ScGeom_ViscElPhys_Basic()],label = 'interactionLoop'),
# Integrate the equation and calculate the new position/velocities...
NewtonIntegrator(damping=0.2, gravity=gravityVector, label='newtonIntegr'),
]

#Timestep
O.dt=.1e-2
qt.View()
O.run(21500)
O.saveTmp('init_channel.xml.bz2')

In the second, i add the inclination with saving evolution of the average kinetic energy of the grains and their average distance in a file and when i open it i found that the kinitic energy=2 in all the simulation, this is my problem.
#!/usr/bin/python
#Import libraries

#Inclination angle of the channel slope
slope=35

#Evolution of the kinetic energies of the flow over the average distance traveled by the grains.
n=7200 # number of spheres

# interactions
O.engines=O.engines+[NewtonIntegrator(damping=0.2, gravity=Vector3(-1*sin(slope),0.0,-1*cos(slope)), label='newtonIntegr'),
PyRunner(iterPeriod=500,command='AddData()'),#save every 500 periode of simulation
PyRunner(iterPeriod=23607,command='Data()'),# end of simulation
]

""" open file"""
fil=open('inclined_channel_40w_.dat', 'w')
fil.write('#Ec\t x \n')

x_g=0
for b in O.bodies:
x_g=x_g+b.state.pos
fil.write('%d %.3f\n' % (e,x_g/n))

#saving simulation
def Data():

#Timestep
O.dt=.1e-2
O.run(23707)
qt.View()
#utils.waitIfBatch()

 Robert Caulk (rcaulk) said on 2020-01-20: #3

Consider creating a minimal working example (MWE ), for increasing the likelihood of assistance. For example, I do not think saving and reloading the simulation in multiple scripts is contributing to your issue. Thus, those "features", and others (is a periodic box even necessary?) should be removed during your creation of the MWE.

But maybe someone will be willing to save both scripts to their computer, edit the file paths for loading and reloading, and run them in correct sequence. But not me ;)

 Jan Stránský (honzik) said on 2020-01-20: #4

What version of Yade and Python you are using?

> e=kineticEnergy()/n

this is int value in Python2 and a float value in Python3. If you use Python2, maybe the value is changing, but it is rounded to the value of 2..

> O.engines=O.engines+[NewtonIntegrator(...)]

be aware, that you already have NewtonIntegrator in your O.engines, you have two integrators after this command, does not sound OK..

> fil.write('%d %.3f\n' % (e,x_g/n))

I guess it crashes on something like "e is not defined"

> for b in O.bodies:
> ...

I would put this line outside the "for b in O.bodies". It just uselessly add the same data 7200 times..

I have tried your code, but the first part "freezes" (is just too slow) around iter 20900 (ironically it should save the final state at iter 21000)

cheers
Jan

 Jan Stránský (honzik) said on 2020-01-20: #5

>> fil.write('%d %.3f\n' % (e,x_g/n))
> I guess it crashes on something like "e is not defined"

ok, sorry, it does not crash, because e IS defined and is math.e (yade uses "from math import *" at startup), Euler's number, 2.718...
If you save it as int (%d, why saving energy as int??), you get the value 2

cheers
Jan

 hafsa (sebbah.hafsa) said on 2020-01-20: #6

1. i am using yade 2018.02b and python3
2. in my case and in my opinion divide the problem in two codes is the better way to facilitate the study and accelerate simulation
3. i used two NewtonIntegrator because i didn't found another command to tilt the system.
4. for the kinetic energy if i use plot.plot the value change over time but when i saved it in file it's always equal 2 even i use (%.3f).
is there any other method to have values of kinetic energy over time apart plot.addData(i=O.iter,e=kineticEnergy()/n,**O.energy) ? Jan Stránský (honzik) said on 2020-01-20: #7

> 3. i used two NewtonIntegrator because i didn't found another command to tilt the system.

maybe
newtonIntegr.gravity = ...
?

Please open a new question if you want further discuss this topic

4. ... but when i saved it [kinetic energy] ...

once more, you are NOT saving kinetic energy, but the value of "e", which is equal to math.e (because "e" is not assign any other value anywhere in your script).

> it's always equal 2 even i use (%.3f).

I used
fil.write('%.3f %.3f\n' % (e,x_g/n))
and the value saved in the file was
2.718

> is there any other method to have values of kinetic energy over time apart plot.addData(i=O.iter,e=kineticEnergy()/n,**O.energy) ?

yes, infinitely many methods. E.g. the one you use (just save what you want to save)

cheers
Jan

PS: ok, the solution:
###
x_g=0
for b in O.bodies:
x_g=x_g+b.state.pos
e = kineticEnergy()/n # !! this was missing in your original script
plot.addData(i=O.iter,e=e,**O.energy) # !! moved outside "for b in O.bodies:"
fil.write('%.3f %.3f\n' % (e,x_g/n))
###

 hafsa (sebbah.hafsa) said on 2020-01-25: #8

Thanks Jan Stránský, that solved my question.

 hafsa (sebbah.hafsa) said on 2020-01-25: #9

thank's Jan,
i used newtonIntegr.gravity=(1*sin(17),0,-1*cos(17)) #to change the gravity