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

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

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()

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

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

>> 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

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) ?

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

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