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

Asked by hafsa on 2020-01-17

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
from yade import qt ,plot

# load the initial congiguration
O.load('/home/hafsa/doctorat d1/init_channel_40w.yade.gz')

#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')

def AddData():
    x_g=0
    for b in O.bodies:
 if isinstance(b.shape,Sphere) and b.shape.radius==0.5:
                plot.addData(i=O.iter,e=kineticEnergy()/n,**O.energy) #energy kenetic
         x_g=x_g+b.state.pos[0]
    fil.write('%d %.3f\n' % (e,x_g/n))

#saving simulation
def Data():
    O.save('/home/hafsa/doctorat d1/inclined_channel_40w.yade.gz')

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

Question information

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

Hello,

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

cheers
Jan

[1] https://www.yade-dem.org/wiki/Howtoask

hafsa (sebbah.hafsa) said : #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
r_Prt=0.5 #radius of the particles
r_Prt_base=0.05 #radius of the particles which made the bunpy base

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):
      O.bodies.append([utils.sphere(center=(i*0.5,j*0.5,1),radius=r_Prt_base,fixed=True,color=(1,1,0))])
      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'),
        PyRunner(iterPeriod=21000,command='AddData()')
 ]

def AddData():
   O.save('/home/hafsa/doctorat d1/init_channel_40w.yade.gz')

#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
from yade import qt ,plot

# load the initial congiguration
O.load('/home/hafsa/doctorat d1/init_channel_40w.yade.gz')

#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')

def AddData():
    x_g=0
    for b in O.bodies:
 if isinstance(b.shape,Sphere) and b.shape.radius==0.5:
                plot.addData(i=O.iter,e=kineticEnergy()/n,**O.energy) #energy kenetic
         x_g=x_g+b.state.pos[0]
    fil.write('%d %.3f\n' % (e,x_g/n))

#saving simulation
def Data():
    O.save('/home/hafsa/doctorat d1/inclined_channel_40w.yade.gz')

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

Robert Caulk (rcaulk) said : #3

Consider creating a minimal working example (MWE [1]), 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 ;)

[1] https://www.yade-dem.org/wiki/Howtoask

Jan Stránský (honzik) said : #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:
> ...
> plot.addData(i=O.iter,e=kineticEnergy()/n,**O.energy) #energy kenetic

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

Please follow Robert's suggestion about MWE.
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 : #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 : #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) ?

Best Jan Stránský (honzik) said : #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:
###
def AddData():
   x_g=0
   for b in O.bodies:
      if isinstance(b.shape,Sphere) and b.shape.radius==0.5:
         x_g=x_g+b.state.pos[0]
   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 : #8

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

hafsa (sebbah.hafsa) said : #9

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