Simulation blocked after O.run in a conditional translation motion

Asked by Rioual on 2020-01-23

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

 TransEngload2 = TranslationEngine(ids=fctIdsWS,translationAxis=[0,-1,0],velocity=10)

 O.engines=O.engines+[TransEngload2]

 while h_WS > hmaxSpheres:

  print 'test01'

  TransEngload2.dead = False

  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

  TransEngload2.dead = True

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

Thank you for your Help,

Best

Vincent

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-29
Last query:
2020-01-29
Last reply:
2020-01-27
Jan Stránský (honzik) said : #1

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

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

Jan Stránský (honzik) said : #3

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

Please read [1] once more.
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 : #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 : #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() ?

Thanks for your reply,

Best Vincent,

Jan Stránský (honzik) said : #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 : #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 : #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 : #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:

Welcome to Yade 2016.06a
TCP python prompt on localhost:9004, auth cookie `skduye'
XMLRPC info provider on http://localhost:21004
Running script working-example.py
End building packing
/home/irstea/myYade/install/lib/x86_64-linux-gnu/yade-2016.06a/py/yade/utils.py:351: UserWarning: Function utils.facetCylinder is deprecated, use geom.facetCylinder instead.
  _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 : #11

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

Jan Stránský (honzik) said : #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)
# 2. readParamsFromTable which reads respective line from the parameter file
# 3. the simulation muse be run using yade-batch, not yade
#
# $ yade-batch --job-threads=1 03-oedometric-test.table03-oedometric-test.py
#

# load parameters from file if run in batch
# default values are used if not run from batch
readParamsFromTable(rMean=.05,rRelFuzz=.3,maxLoad=1e6,minLoad=1e4)
# make rMean, rRelFuzz, maxLoad accessible directly as variables later
from yade.params.table import *

# create box with free top, and ceate loose packing inside the box
from yade import pack, plot
O.bodies.append(geom.facetBox((.5,.5,.5),(.5,.5,.5),wallMask=31))
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 following checkUnbalanced, unloadPlate and stopUnloading functions are all called by the 'checker'
# (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
 O.bodies.append(wall(max([b.state.pos[2]+b.shape.radius for b in O.bodies
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
 O.engines=O.engines+[PyRunner(command='addPlotData()',iterPeriod=200)]
 # next time, do not call this function anymore, but the next one (unloadPlate) instead
 checker.command='unloadPlate()'

def unloadPlate():
 # if the force on plate exceeds maximum load, start unloading
 if abs(O.forces.f(plate.id)[2])>maxLoad:
  plate.state.vel*=-1
  # next time, do not call this function anymore, but the next one (stopUnloading) instead
  checker.command='stopUnloading()'

def stopUnloading():
 if abs(O.forces.f(plate.id)[2])<minLoad:
  # 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()

def addPlotData():
 if not isinstance(O.bodies[-1].shape,Wall):
  plot.addData(); return
 Fz=O.forces.f(plate.id)[2]
#
 plot.addData(Fz=Fz,w=plate.state.pos[2]-plate.state.refPos[2],unbalanced=unbalancedForce(),i=O.iter)

# 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
from yade import qt
from yade import ymport
from yade.gridpfacet import *

###########################################
## PhysicalParameters
Density=2400
frictionAngle=radians(35)
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))
from yade import ymport
#******************************************
global TransEngload2
#******************************************
# Importation de la roue

global fctIdsWS

fctIdsWS=O.bodies.append( facetCylinder(center=(0.5,0.5,1.2), radius=0.3,
height=0.2, orientation=Quaternion((1, 0, 0), 0), segmentsNumber=10,
wallMask=7, angleRange=None, closeGap=False, radiusTopInner=-1,
radiusBottomInner=-1, material=facetMat) )

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

from yade import qt
qt.View()

def kinematics_WS():

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

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

 TransEngload2 = TranslationEngine(ids=fctIdsWS,translationAxis=[0,0,-1],velocity=10)

 O.engines=O.engines+[TransEngload2]

 while h_WS > hmaxSpheres:

  print 'test01'

  TransEngload2.dead = False

  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

  TransEngload2.dead = True

 else:

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

  while (amTOT< 1e-10):

   O.engines=O.engines+[PyRunner(command='addTorque()',iterPeriod=1)]

   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
# else: return (hMax+O.bodies[idHMax].shape.radius)

 hmaxSpheres = hMax+O.bodies[idHMax].shape.radius

 return h_WS,hmaxSpheres

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

def addTorque():

 increment_T= 0.001

 for Idfacets in fctIdsWS:
  O.forces.addT(Idfacets,(0,increment_T,0))
  imposed_T = O.forces.permT(Idfacets)
 return imposed_T

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

Jan Stránský (honzik) said : #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 : #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!

Thanks for your feedback,

Best

V.

Best Jan Stránský (honzik) said : #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 : #18

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