# Simulation blocked after O.run in a conditional translation motion

Hello,

This mistake may be obvious but I can't get it.
I am studying the translation motion of a boundary WS (using TranslationEngine) as long as an conditional statement is fulfilled
(the height of the boundary h_WS must be smaller than the maximum height of the packing hmaxSpheres).
I use O.run(1,True) to let the engine work for one time step but the simulation seems to be frozen at the "O.run(1,True)" statement.

****************************************************************************************************
[...
PyRunner(command='kinematics_WS()',realPeriod=1,label='kine'),
...
]

def kinematics_WS():

h_WS = calc_h()[0]
hmaxSpheres = calc_h()[1]

print 'test0','h_WS=',h_WS,'hmaxSpheres=',hmaxSpheres

while h_WS > hmaxSpheres:

print 'test01'

print 'test02'

O.run(1,True) ***************frozen here!***************************

# O.engines=O.engines+[PyRunner(command='calc_h()')]

print 'test03'

h_WS = calc_h()[0]
hmaxSpheres = calc_h()[1]

print 'test1','h_WS=',h_WS,'hmaxSpheres=',hmaxSpheres

else:
.....

****************************************************************************************************************
Output:
test0 h_WS= 130.0 hmaxSpheres= 100.369406791
test01
test02
.......frozen
The boundary does not seem to move even if i put a higher step number in O.run()
*****************************************************************************************************************

Best

Vincent

## Question information

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

Hello,
please provide a complete script [1], otherwise the best we can do is guessing..
thanks
Jan

 Jan Stránský (honzik) said on 2020-01-24: #3

> fctIdsWS=O.bodies.append(ymport.stl('WS-FR-1.stl',color=(1,0,0),material=facetMat))

By complete script I meant a code which we can try - which is not possible in this case..

So another information request :-) freezing at the very first step can just mean that Yade is computing.
For how long have you let it be frozen? a few seconds? one hour? a day?
How many particles (spheres, facets, ... ) do you have in your simulation? (it may be that it is just too many particles to be computed quickly)

cheers
Jan

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

Hello,

in this specific case, you have O.run() in the script and then O.run(1,True).. (calling run on already running simulation)
If I add O.run();O.wait() in the initial phase, the code will run quickly, not reaching the O.run(1,True) branch..

cheers
Jan

 Rioual (francois-rioual-v) said on 2020-01-24: #6

Hello,

But I have a O.pause() in the function stopUnloading() in order to stop the process of building the packing (first part).
So this should desactivate the first O.run(); No ?
I think I am confused in dealing with O.run() ?

Best Vincent,

 Jan Stránský (honzik) said on 2020-01-24: #7

O.pause() would pause the running, yes, but before it is reached, python continues to execute other lines one by one. Therefore you need O.wait() to force python wait for the O.pause()
Jan

 Rioual (francois-rioual-v) said on 2020-01-24: #8

Hello Jan,

For me it does not change anything adding O.wait().
It is still blocking at the same line !

regards

V.

 Jan Stránský (honzik) said on 2020-01-24: #9

strange, because when I use O.run();O.wait() instead just the first O.run(), the O.run(1,True) code is not executed at all..
what version of Yade do you use?

 Rioual (francois-rioual-v) said on 2020-01-24: #10

......Yes, it is still blocking. My yade version is Yade 2016.06a apparently.

I added a print before the pause in the script:
print 'End building packing'
O.pause()

**************************************************************************
Here is my output:

TCP python prompt on localhost:9004, auth cookie `skduye'
XMLRPC info provider on http://localhost:21004
Running script working-example.py
End building packing
_deprecatedUtilsFunction('facetCylinder','geom.facetCylinder')
[[ ^L clears screen, ^U kills line. F12 controller, F11 3d view (use h-key for showing help), F10 both, F9 generator, F8 plot. ]]
[0;34mYade [[1;34m1[0;34m]: [0mWARN /home/irstea/myYade/trunk/trunk-2016.06a/pkg/common/GravityEngines.cpp:19 action: GravityEngine is deprecated, consider using Newton::gravity instead (unless gravitational energy has to be tracked - not implemented with the newton attribute).
test0 h_WS= 1.3 hmaxSpheres= 0.443123580808
test01
test02
....frozen....
************************************************************************

Regards

Vincent

 Jan Stránský (honzik) said on 2020-01-24: #11

To be sure, please post once more exactly the code you are currently using
thanks
Jan

 Jan Stránský (honzik) said on 2020-01-25: #13

I have tried once more the code using Yade 2018.02b on Ubuntu 18.04 with the same result - the problematic code is not reached at all..

Ibelow is the exact code I used (I had some problem with formatting the copy-pasted code and got syntax errors)

If the code is problematic for you, then the yade version contributes to the problem (and personally I don't want to search for problems with 2016 version..)

cheers
Jan

#####
# gravity deposition, continuing with oedometric test after stabilization
# shows also how to run parametric studies with yade-batch

# The components of the batch are:
# 1. table with parameters, one set of parameters per line (ccc.table)
#
#

# load parameters from file if run in batch
# default values are used if not run from batch
# make rMean, rRelFuzz, maxLoad accessible directly as variables later

# create box with free top, and ceate loose packing inside the box
sp=pack.SpherePack()
sp.makeCloud((0,0,0),(1,1,1),rMean=rMean,rRelFuzz=rRelFuzz)
sp.toSimulation()

O.engines=[
ForceResetter(),
# sphere, facet, wall
InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb(),Bo1_Wall_Aabb()]),
InteractionLoop(
# the loading plate is a wall, we need to handle sphere+sphere, sphere+facet, sphere+wall
[Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom(),Ig2_Wall_Sphere_ScGeom()],
[Ip2_FrictMat_FrictMat_FrictPhys()],
[Law2_ScGeom_FrictPhys_CundallStrack()]
),
NewtonIntegrator(gravity=(0,0,-9.81),damping=0.5),
# the label creates an automatic variable referring to this engine
# we use it below to change its attributes from the functions called
PyRunner(command='checkUnbalanced()',realPeriod=2,label='checker'),
]
O.dt=.5*PWaveTimeStep()

# (the last engine) one after another; this sequence defines progression of different stages of the
# simulation, as each of the functions, when the condition is satisfied, updates 'checker' to call
# the next function when it is run from within the simulation next time

# check whether the gravity deposition has already finished
# if so, add wall on the top of the packing and start the oedometric test
def checkUnbalanced():
# at the very start, unbalanced force can be low as there is only few contacts, but it does not mean the packing is stable
if O.iter<5000: return
# the rest will be run only if unbalanced is < .1 (stabilized packing)
if unbalancedForce()>.1: return
# add plate at the position on the top of the packing
# the maximum finds the z-coordinate of the top of the topmost particle
if isinstance(b.shape,Sphere)]),axis=2,sense=-1))
global plate # without this line, the plate variable would only exist inside this function
plate=O.bodies[-1] # the last particles is the plate
# Wall objects are "fixed" by default, i.e. not subject to forces
# prescribing a velocity will therefore make it move at constant velocity (downwards)
plate.state.vel=(0,0,-.1)
# start plotting the data now, it was not interesting before
# next time, do not call this function anymore, but the next one (unloadPlate) instead

plate.state.vel*=-1

# O.tags can be used to retrieve unique identifiers of the simulation
# if running in batch, subsequent simulation would overwrite each other's output files otherwise
# d (or description) is simulation description (composed of parameter values)
# while the id is composed of time and process number
# plot.saveDataTxt(O.tags['d.id']+'.txt')
print 'End building packing'
O.pause()

if not isinstance(O.bodies[-1].shape,Wall):
Fz=O.forces.f(plate.id)[2]
#

# besides unbalanced force evolution, also plot the displacement-force diagram
#plot.plots={'i':('unbalanced',),'w':('Fz',)}
#plot.plot()

O.run();O.wait()
#*********************************************************************************************************************
# Let the object penetrate the packing

#Start penetration

import time

###########################################
## PhysicalParameters
Density=2400
tc = 0.001
en = 0.05
et = 0.05

## Import wall's geometry
facetMat=O.materials.append(ViscElMat(frictionAngle=frictionAngle,tc=tc,
en=en, et=et)) # **params sets kn, cn, ks, cs
sphereMat=O.materials.append(ViscElMat(density=Density,frictionAngle=frictionAngle,tc=tc,en=en,et=et))
#******************************************
#******************************************
# Importation de la roue

global fctIdsWS

height=0.2, orientation=Quaternion((1, 0, 0), 0), segmentsNumber=10,

## Timestep
O.dt=.2*tc

## Engines
O.engines=[
## Resets forces and momenta the act on bodies
ForceResetter(),
## Using bounding boxes find possible body collisions.
#InsertionSortCollider(),
InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb(),Bo1_GridConnection_Aabb()]),
## Interactions
InteractionLoop(
# the loading plate is a facet, we need to handle sphere+sphere, sphere+facet
[Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom(),Ig2_Wall_Sphere_ScGeom(),Ig2_Sphere_GridConnection_ScGridCoGeom()],
[Ip2_FrictMat_FrictMat_FrictPhys()],
[Law2_ScGeom_FrictPhys_CundallStrack()]
),
## Apply gravity
GravityEngine(gravity=[0,0,-9.81]),
## Cundall damping must been disabled!
NewtonIntegrator(damping=0),

## Apply kinematics to wheel
PyRunner(command='kinematics_WS()',realPeriod=1,label='kine'),
]

qt.View()

def kinematics_WS():

h_WS = calc_h()[0]
hmaxSpheres = calc_h()[1]

print 'test0','h_WS=',h_WS,'hmaxSpheres=',hmaxSpheres

while h_WS > hmaxSpheres:

print 'test01'

print 'test02'

######################### Here the simulation seems to be blocked#################################################################

O.run(1,True)

# O.engines=O.engines+[PyRunner(command='calc_h()')]

print 'test03'

h_WS = calc_h()[0]
hmaxSpheres = calc_h()[1]

print 'test1','h_WS=',h_WS,'hmaxSpheres=',hmaxSpheres

else:

amTOT = sum((O.bodies[facetid].state.angMom)[1] for facetid in fctIdsWS)

while (amTOT< 1e-10):

amTOT = sum((O.bodies[facetid].state.angMom)[1] for facetid in fctIdsWS)

else:

# Stop simulation et measurement of the shear torque and cohesion C_WS = 3* (imposed_T)/(2*Pi*(R0^3-R1^3))
print 'C_WS=',C_WS

print '********End of the calculation for Cohesion***********' #
O.pause()

#*********************************************************************************************************************************
#This function calculate the height of the boundary h_WS and the maximum height of the packing hmax_Spheres
def calc_h():

###########################################
#Calculate h_WS

minZ = 1e99
maxZ = -1e99

for facet in fctIdsWS:
# print 'facet=', facet
facet= O.bodies[facet]
vs = [facet.state.pos + facet.state.ori * v for v in facet.shape.vertices] # vertices in global coord system
# print 'vs=',vs

minZ = min(minZ,min(v[2] for v in vs))
maxZ = max(maxZ,max(v[2] for v in vs)) #### print 'maxZ=',maxZ,'minZ=',minZ

h_WS = maxZ

###########################################
#Calculate hmax_Spheres

idHMax=0 # on definit une variable pour identifier quel corps a la hauteur la plus haute
hMax=0.0 # initialisation de la hauteur a zero
for i in O.bodies: # on parcours tout les corps du systeme
h=i.state.pos[2] # on extrait la position selon l'axe y
if (type(i.shape).__name__ == 'Sphere' and h>hMax): #si le solide est une sphere et sa position est plus haut que hmax
hMax=h # on enregistre sa hauteur
idHMax=i.id # on garde en memoire de quel corps il s'agit
# if (idHMax == 0): return rMax # valeur de retour par default, rMAX

return h_WS,hmaxSpheres

#*****************************************************************************************************************************
#*****************************************************************************************************************************
# Add incremental torque on the WS boundary

increment_T= 0.001

for Idfacets in fctIdsWS:
imposed_T = O.forces.permT(Idfacets)
return imposed_T

#*****************************************************************************************************************************
#####

 Jan Stránský (honzik) said on 2020-01-27: #15

Sorry, my bad, I had to do something wrong, I also get frozen simulation..

probably it is because there is O.run inside a function, which is called from PyRunner (it is there from the OP, sorry for not noticing it..)
Could you just call the function "externally", not from O.engines?

Jan

 Rioual (francois-rioual-v) said on 2020-01-27: #16

Hello Jan,

The function 'kinematic_WS' is ruling all the different dynamic operations I have to do with the wheel in
this code: A translation of the wheel under a certain condition then a rotation of the wheel under a certain
condition.
I am not a yade expert and I really do not know how to express that differently and to call the function
externally!

Best

V.

 Jan Stránský (honzik) said on 2020-01-27: #17

"extenrally" I meant:
- delete the PyRunner(command='kinematics_WS()',...) from O.engines
- call kinematics_WS() somewhere around the end of the script.

simply, if you need O.step or O.run in a function, do not use it in a PyRunner

Jan

 Rioual (francois-rioual-v) said on 2020-01-29: #18

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