Issues in using psd to create a cloud of spheres

Hi all,

I have a code that generate a cloud of spheres and compress the cloud to a certain porosity. It works well when I make the cloud with specific radius and rRelFuzz. However, when I try to make the cloud with specific size distribution using psdSizes and psdCumm, the code run but the spheres will freeze and not compact. Can you please help and let me know what went wrong? My code is below.

Thank you,
Othman

----------

from yade import pack

targetp = .3 ##specify the targeted porosity

##specimen geometry
radiuscyl=.05
heightcyl=0.3

# material parameters

O.materials.append(FrictMat(young = 5e10, poisson = 0.15,frictionAngle = atan(.2), density=1920))

############################ spheres #############################
sp=pack.SpherePack()
psdSizes=[0,.00238,.00475,.00952,.0127]
psdCumm=[0,0.0173,.0289,.9469,1]
sp.makeCloud((0,0,0),(.3,.3,1.5),psdSizes=psdSizes,psdCumm=psdCumm)
#sp.makeCloud((0,0,0),(.3,.3,1.5),rMean=.00475) #if I uncomment this line and comment the above 3 lines, the spheres will not move
#### cylinder extraction
pred=pack.inCylinder((.2,.2,0),(.2,.2,heightcyl),radiuscyl)
spFilter=filterSpherePack(pred,sp,Material=Material, returnSpherePack=True)
spFilter.toSimulation()
############################ facets #############################

facets=geom.facetCylinder((.2,.2,heightcyl/2),radiuscyl,heightcyl,segmentsNumber=50,wallMask=4)

cylinder=O.bodies.append(facets)
yade.qt.View()

##creating disks

d1=geom.facetCylinder((.2,.2,heightcyl),radiuscyl-.002,0,segmentsNumber=50,wallMask=1)
d2=geom.facetCylinder((.2,.2,0),radiuscyl-.002,0,segmentsNumber=50,wallMask=1)

disk1IDs= O.bodies.append(d1)
disk2IDs= O.bodies.append(d2)

for i in disk1IDs:
 body= O.bodies[i]
 body.state.vel = (0,0,-10)

for n in disk2IDs:
 body= O.bodies[n]
 body.state.vel = (0,0,10)

############################ compaction #############################
O.dt=.5*utils.PWaveTimeStep()
enlargeFactor=1.5
O.engines=[
 ForceResetter(),
 InsertionSortCollider([
  Bo1_Sphere_Aabb(aabbEnlargeFactor=enlargeFactor,label='bo1s'),
  Bo1_Facet_Aabb()
 ]),
 InteractionLoop(
  [
   Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=enlargeFactor,label='ss2d3dg'),
   Ig2_Facet_Sphere_ScGeom(),
  ],
  [

   Ip2_FrictMat_FrictMat_FrictPhys(),
   Ip2_FrictMat_FrictMat_FrictPhys(),
  ],
  [

   Law2_ScGeom_FrictPhys_CundallStrack(),
  ],
 ),

 NewtonIntegrator(damping=.3),

 PyRunner(iterPeriod=500,command='stop()'),

]
O.run()
# reset interaction detection enlargement
bo1s.aabbEnlargeFactor=ss2d3dg.interactionDetectionFactor=1.0

def stop():
   if utils.porosity()<targetp:
 O.pause()
 print 'Finished'
 for i in disk1IDs: O.bodies.erase(i)
 for i in disk2IDs: O.bodies.erase(i)
 for i in cylinder: O.bodies.erase(i)
 print ('unbalanced forces = ', utils.unbalancedForce())
 print ('cylinder dimensions = ', utils.aabbDim())
 print ('porosity = ', utils.porosity())
 print ('No. of spheres = ', len (O.bodies))