gravity deposition of clumps

if there is a better way to generate the geometry I would be also interested to hear it

Hello,

concerning plot "error":

> but end up with error
> UserWarning: No labelled objects found. Use label='...' kwarg on individual plots.

it is not an error, just a warning. The problem is probably that at the plot creation time, there are no plot data (but plot.plots expects some).

using O.run(1000,True) before plot.plot() command solves it

Jan

Thank you very much Jan for the quick reply! Indeed O.run was needed. Unfortunately the gravity doesnt seem to influence the clumps.

Ive update the code

######################################################################
# A script for creating a dense packing of clumps
#
# Each aggragate is a dense packing, but macroscopically the packing
# is loose
######################################################################
def getClumpInfo():
 c = []
 for b in O.bodies:
  #print(b.isClump)
  if b.isClump:
   c.append(b)
   #print(b.shape)
   print('Clump ',b.id,' has following members:')
   keys = b.shape.members.keys()
   for ii in range(0,len(keys)):
    print('- Body ',keys[ii])
 #print(c)
 return c

def setClumpInfo(disp=False):
 c = []
 aggNum = -1
 for b in O.bodies:
  #print(b.isClump)

  if b.isClump:
   aggNum+=1
   c.append(b)
   #print(b.shape)
   if disp: print('Clump ',b.id,' has following members:')
   keys = b.shape.members.keys()
   for ii in range(0,len(keys)):
    if disp: print('- Body ',keys[ii])
    #O.bodies[keys[ii]].agglomerate = b.id # tell each particle who is its agglomerate
    O.bodies[keys[ii]].agglomerate = aggNum # tell each particle who is its agglomerate
 print(str(len(c))+ ' ' + str(aggNum))
 return len(c) #number of clumps

# if the unbalanced forces goes below .05, the packing
# is considered stabilized, therefore we stop collected
# data history and stop
def checkUnbalanced():
 from yade import pack, plot
 if unbalancedForce()<.05:
  O.pause()
  plot.saveDataTxt('bbb.txt.bz2')
  # plot.saveGnuplot('bbb') is also possible

# collect history of data which will be plotted
def addPlotData():
 from yade import pack, plot
 # each item is given a names, by which it can be the unsed in plot.plots
 # the **O.energy converts dictionary-like O.energy to plot.addData arguments
 plot.addData(i=O.iter,unbalanced=unbalancedForce(),**O.energy)

def gravityDeposition():
 # gravity deposition in box, showing how to plot and save history of data,
 # and how to control the simulation while it is running by calling
 # python functions from within the simulation loop

 # import yade modules that we will use below
 from yade import pack, plot

 # create rectangular box from facets
 O.bodies.append(geom.facetBox(center,extend,wallMask=31))

 O.engines=[
    ForceResetter(),
    InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()]),
    InteractionLoop(
    # handle sphere+sphere and facet+sphere collisions
    [Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom()],
    [Ip2_FrictMat_FrictMat_FrictPhys()],
    [Law2_ScGeom_FrictPhys_CundallStrack()]
    ),
    NewtonIntegrator(gravity=(0,0,-9.81),damping=0.4),
    # call the checkUnbalanced function (defined below) every 2 seconds
    PyRunner(command='checkUnbalanced()',realPeriod=2),
    # call the addPlotData function every 200 steps
    PyRunner(command='addPlotData()',iterPeriod=100)
 ]
 O.dt=.5*PWaveTimeStep()

 # enable energy tracking; any simulation parts supporting it
 # can create and update arbitrary energy types, which can be
 # accessed as O.energy['energyName'] subsequently
 O.trackEnergy=True

 #O.run(500,True)

 raw_input('test')

 # define how to plot data: 'i' (step number) on the x-axis, unbalanced force
 # on the left y-axis, all energies on the right y-axis
 # (O.energy.keys is function which will be called to get all defined energies)
 # None separates left and right y-axis
 plot.plots={'i':('unbalanced',None,O.energy.keys)}

 O.run(1000,True)

 # show the plot on the screen, and update while the simulation runs
 plot.plot()

from yade import pack,export,ymport
import random
random.seed(1) # to make colors always the same

numClumps = 2

minCorner = (0,0,0)
maxCorner = (7e-7,7e-7,7e-7)

#O.periodic=True
sp = pack.SpherePack()

# load clumps
ymport.textClumps('/tmp/clump0.txt')

c0 = pack.SpherePack()
c0.fromSimulation()

ymport.textClumps('/tmp/clump1.txt')
#c1 = getClumpInfo()
c1 = pack.SpherePack()
c1.fromSimulation()

O.switchScene(); O.resetThisScene() #####!!!!!!!

print('Creating Clump Cloud')
test = sp.makeClumpCloud(minCorner, maxCorner, [c0,c1], periodic=True, num = 1000)#, periodic=True, num=-1, seed=1)
print(test)

sp.toSimulation()

dim=utils.aabbExtrema()
dim=utils.aabbExtrema()
xinf=dim[0][0]
xsup=dim[1][0]
xdim = xsup-xinf
X=xinf+(xdim)/2.
yinf=dim[0][1]
ysup=dim[1][1]
ydim = ysup-yinf
Y=yinf+(ydim)/2.
zinf=dim[0][2]
zsup=dim[1][2]
zdim = zsup-zinf
Z=zinf+(zdim)/2.

center = Vector3(X,Y,Z) #center of the packing
extend = Vector3(xdim,ydim,zdim) #extends of the packing

raw_input('Before gravity')
gravityDeposition()
raw_input('After gravity')
quit()

setClumpInfo()
#print(enumerate(sp))

# save the result, including information of agglomerates which the particle belongs to
export.textExt('/tmp/divided.txt','x_y_z_r_attrs',attrs=['b.agglomerate'])

print('Number of particles = ' , len(O.bodies))

try:
 from yade import qt
 qt.View()
except:
 pass

Note that for https://yade-dem.org/doc/tutorial-examples.html#gravity-deposition (and probably your example as well ?), you also have
"
yade.plot: creating fake plot, since there are no y-data yet
"

(hopefully clear enough ?) before the warning "UserWarning: No labelled objects found. Use label='...' kwarg on individual plots" you mentioned.

Hi Jérôme, thank you for the reply.

Since adding O.run the 'No labelled objects found' error doesnt appear.

However I am still having issues where the gravity I set doesnt influence the clumps

please check that the code you provide is a working code. Trying your updated version I got
NameError: global name 'extend' is not defined
thanks
Jan

Hi Jan,
It should work as extend was a global variable.
In case it still doesnt work, Ive moved the section where extend is used to be a part of gravityDeposition()

######################################################################
# A script for creating a dense packing of clumps
#
# Each aggragate is a dense packing, but macroscopically the packing
# is loose
######################################################################
def getClumpInfo():
 c = []
 for b in O.bodies:
  #print(b.isClump)
  if b.isClump:
   c.append(b)
   #print(b.shape)
   print('Clump ',b.id,' has following members:')
   keys = b.shape.members.keys()
   for ii in range(0,len(keys)):
    print('- Body ',keys[ii])
 #print(c)
 return c

def setClumpInfo(disp=False):
 c = []
 aggNum = -1
 for b in O.bodies:
  #print(b.isClump)

  if b.isClump:
   aggNum+=1
   c.append(b)
   #print(b.shape)
   if disp: print('Clump ',b.id,' has following members:')
   keys = b.shape.members.keys()
   for ii in range(0,len(keys)):
    if disp: print('- Body ',keys[ii])
    #O.bodies[keys[ii]].agglomerate = b.id # tell each particle who is its agglomerate
    O.bodies[keys[ii]].agglomerate = aggNum # tell each particle who is its agglomerate
 print(str(len(c))+ ' ' + str(aggNum))
 return len(c) #number of clumps

# if the unbalanced forces goes below .05, the packing
# is considered stabilized, therefore we stop collected
# data history and stop
def checkUnbalanced():
 from yade import pack, plot
 if unbalancedForce()<.05:
  O.pause()
  plot.saveDataTxt('bbb.txt.bz2')
  # plot.saveGnuplot('bbb') is also possible

# collect history of data which will be plotted
def addPlotData():
 from yade import pack, plot
 # each item is given a names, by which it can be the unsed in plot.plots
 # the **O.energy converts dictionary-like O.energy to plot.addData arguments
 plot.addData(i=O.iter,unbalanced=unbalancedForce(),**O.energy)

def gravityDeposition():
 # gravity deposition in box, showing how to plot and save history of data,
 # and how to control the simulation while it is running by calling
 # python functions from within the simulation loop

 # import yade modules that we will use below
 from yade import pack, plot

 dim=utils.aabbExtrema()
 dim=utils.aabbExtrema()
 xinf=dim[0][0]
 xsup=dim[1][0]
 xdim = xsup-xinf
 X=xinf+(xdim)/2.
 yinf=dim[0][1]
 ysup=dim[1][1]
 ydim = ysup-yinf
 Y=yinf+(ydim)/2.
 zinf=dim[0][2]
 zsup=dim[1][2]
 zdim = zsup-zinf
 Z=zinf+(zdim)/2.

 center = Vector3(X,Y,Z) #center of the packing
 extend = Vector3(xdim,ydim,zdim) #extends of the packing

 # create rectangular box from facets
 O.bodies.append(geom.facetBox(center,extend,wallMask=31))

 O.engines=[
    ForceResetter(),
    InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()]),
    InteractionLoop(
    # handle sphere+sphere and facet+sphere collisions
    [Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom()],
    [Ip2_FrictMat_FrictMat_FrictPhys()],
    [Law2_ScGeom_FrictPhys_CundallStrack()]
    ),
    NewtonIntegrator(gravity=(0,0,-9.81),damping=0.4),
    # call the checkUnbalanced function (defined below) every 2 seconds
    PyRunner(command='checkUnbalanced()',realPeriod=2),
    # call the addPlotData function every 200 steps
    PyRunner(command='addPlotData()',iterPeriod=100)
 ]
 O.dt=.5*PWaveTimeStep()

 # enable energy tracking; any simulation parts supporting it
 # can create and update arbitrary energy types, which can be
 # accessed as O.energy['energyName'] subsequently
 O.trackEnergy=True

 #O.run(500,True)

 raw_input('test')

 # define how to plot data: 'i' (step number) on the x-axis, unbalanced force
 # on the left y-axis, all energies on the right y-axis
 # (O.energy.keys is function which will be called to get all defined energies)
 # None separates left and right y-axis
 plot.plots={'i':('unbalanced',None,O.energy.keys)}

 O.run(1000,True)

 # show the plot on the screen, and update while the simulation runs
 plot.plot()

from yade import pack,export,ymport
import random
random.seed(1) # to make colors always the same

numClumps = 2

minCorner = (0,0,0)
maxCorner = (7e-7,7e-7,7e-7)

#O.periodic=True
sp = pack.SpherePack()

# load clumps
ymport.textClumps('/tmp/clump0.txt')

c0 = pack.SpherePack()
c0.fromSimulation()

ymport.textClumps('/tmp/clump1.txt')
#c1 = getClumpInfo()
c1 = pack.SpherePack()
c1.fromSimulation()

O.switchScene(); O.resetThisScene() #####!!!!!!!

print('Creating Clump Cloud')
test = sp.makeClumpCloud(minCorner, maxCorner, [c0,c1], periodic=True, num = 1000)#, periodic=True, num=-1, seed=1)
print(test)

sp.toSimulation()

raw_input('Before gravity')
gravityDeposition()
raw_input('After gravity')
quit()

setClumpInfo()
#print(enumerate(sp))

# save the result, including information of agglomerates which the particle belongs to
export.textExt('/tmp/divided.txt','x_y_z_r_attrs',attrs=['b.agglomerate'])

print('Number of particles = ' , len(O.bodies))

try:
 from yade import qt
 qt.View()
except:
 pass

sorry, I just wrongly copied your code..
anyway I got
...
ymport.textClumps('/tmp/clump0.txt')
...
IOError: [Errno 2] No such file or directory: '/tmp/clump0.txt'

Jan

Sorry!
The two files clump0.txt and clump1.txt I have pasted here
https://pastebin.com/APtLQdvH

I found another thread which used makeclumpcloud and gravity (https://answers.launchpad.net/yade/+question/290933)

I added the following two lines to see if they would help
O.materials.append(FrictMat(young=1e10, poisson=0.28, frictionAngle=0.2, density=2600.))
O.bodies.updateClumpProperties(discretization=10)# correct mass, volume, inertia!!

however the clumps still dont fall

I've also pasted clump0.txt here https://pastebin.com/WY703ZLs

and clump1.txt here: https://pastebin.com/KNacYF9N

> however the clumps still dont fall

please provide us with the info how you know this (visually from GUI, printing displacement ...)

I have tried

clumps = [b for b in O.bodies if b.isClump]
for c in clumps: print c.displ()

with non-zero result, so the clumps clearly DO fall..

cheers
Jan

Hi Jan,
Im looking at the gui, in the 3D viewport the clumps dont move, I left it running for 10 minutes or so.
If the clumps are moving I guess the timestep is too small?
many thanks
Jesse

> Im looking at the gui

this might be misleading, always try to get a numeric evidence for your conclusions..

> I guess the timestep is too small?

might be. Note that for rigid clumps you can use larger timestep than PWaveTimeStep as the mass of a clump is higher than that of a single sphere..

Jan

>you can use larger timestep than PWaveTimeStep as the mass of a clump is higher than that of a single sphere

Note that GlobalStiffness time stepper accounts for that, so it should give automaticaly a larger timestep.
B

Thank you Jan for the reply,
Is there a general guidline as to the maxmium for O.dt?

Ive increased the timestep by 5 times and the 3 or 4 of the clumps now do move, but in opposite directions. Im not sure how to quanitfy this however..

https://yade-dem.org/doc/formulation.html#stability-considerations

Thank you for the reply. Ill try to calculate the maximum stable time step

Unfortunately, regardless of my choice of timestep, the clumps do not settle. Could there be another reason?

> Is there a general guidline as to the maxmium for O.dt?

yes, do not use higher time step than the critical one :-)

> 3 or 4 of the clumps now do move, but in opposite directions.

what does "opposite directions" means?
I have tried it and all clumps fall uniformly in -z direction, as expected..

Jan

Hi Jan,
What did you use for a timestep??
Did you modify the script?
many thanks
Jesse

> Did you modify the script?

Yes, I have modified the script (e.g. I have deleted raw_input which has no place in MWE)

> What did you use for a timestep??

for testing I used very large time step..

> Unfortunately, regardless of my choice of timestep, the clumps do not settle. Could there be another reason?

Please be more specific about "do not settle"

cheers
Jan

> Please be more specific about "do not settle"

I mean the clumps do not fall uniformly in -z direction

Could you specify the timestep you used?

> I mean the clumps do not fall uniformly in -z direction

In my case, I got perfect uniform free fall in -z direction..

> Could you specify the timestep you used?

e.g. O.dt=5000*PWaveTimeStep(), but the value does not influence the uniform free fall

Jan

Hi Jan,
Thank you for the reply. I also tried O.dt=5000*PWaveTimeStep(), and the clumps do fall uniformly in the -z direction. However, they do not stop when reaching the face of the facetBox, they seem just to disappear(see below). The likely cause is the time step is too large?

I lowered then to O.dt=1000*PWaveTimeStep(). The clumps move lowly in the -z direction, however, now when the clumps reach the lower z face of the facetbox, they do not stop there. Instead they move very quickly in a direction (mostly the +z direction (not disappearing then)).

Could this perhaps be an issue with the colliders?
InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()]),
or the interaction loop, or the dampening of the system?

Is there a way to quantize this?

many thanks
Jesse

> they seem just to disappear
> they move very quickly

as you already guessed, the time step is too large (in both cases).
Jan