Imported spheres go through the walls

Dear all,

I find a question that is weird to me. This question is followed by my previous one [1].
I've created and saved a stress-free dense packing with a script (MWE 1) and then I load it with another script (MWE 2) and run isotropic compression.
While I was running the isotropic compression, some of the imported spheres went through the walls. To fix this problem, I've tried to decrease the timestep by setting a lower safety factor (0.8 --> 0.3) and also increase the stiffness of both particles and walls (1e8 --> 1e10). Unfortunately, both methods don't work. I have tried ever higher
So I ran another simulation (MWE 3) to check if this also happens a non-imported loose packing where I generate particles with makeCloud with identical material properties (1e8) and timestep (factor 0.8). The process of isotropic compression is successful as I expected where no particles going through the walls. For this reason, I don't think the particle penetration issue is related to stiffness or timestep.
For the imported dense packing, I imagine it's also a "loose" packing since it's stress-free and there is no overlap between particles.
Under this situation, why would I get particle penetration issue with imported dense packing?
Does someone have any idea about it? I've been struggling with this for a while...
Thanks in advance!

The MWEs are provided below.

[1] https://answers.launchpad.net/yade/+question/694379

MWE 1
##########################################
######## Stress-free packing generation ############
#########################################

from yade import pack,plot,export
import numpy as np
sp=pack.SpherePack()
O.periodic=True
a=5
b=5
c=5
RADIUS_dimension=0.0250
OUT='mypacking'
length=a*(2*RADIUS_dimension)
height=b*(2*RADIUS_dimension)
width=c*(2*RADIUS_dimension)
thickness=RADIUS_dimension

### Particle size distribution
psdSizes=[.0232,.0254,.0276]
psdCumm=[0,0.5,1]

### friction angles
spFRIC=0.5

### boundary conditions
PI=1.e5

### material properties
O.materials.append(CohFrictMat(isCohesive=True,density=2500,young=1e8,poisson=0.5,frictionAngle=spFRIC,normalCohesion=-1,shearCohesion=-1,label='sphereMat'))

O.cell.hSize=Matrix3(length,0,0,0,3*height,0,0,0,width)
sp.makeCloud((0,0,0),(length,3*height,width), psdSizes = psdSizes, psdCumm = psdCumm, periodic=True, seed =1)
O.bodies.append([utils.sphere(s[0],s[1],color=(1,1,0),material='sphereMat') for s in sp])

O.engines=[
 ForceResetter()
 ,InsertionSortCollider([Bo1_Box_Aabb(),Bo1_Sphere_Aabb()],verletDist=-0.1,allowBiggerThanPeriod=True)
 ,InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom6D()],
  [Ip2_CohFrictMat_CohFrictMat_CohFrictPhys()],
  [Law2_ScGeom6D_CohFrictPhys_CohesionMoment()]
 )
 ,GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=1,timestepSafetyCoefficient=0.8,defaultDt=-1)
 ,PeriTriaxController(dynCell=True,mass=10,maxUnbalanced=1e-3,relStressTol=1e-4,stressMask=7,goal=(-PI,-PI,-PI),globUpdate=1,maxStrainRate=(10,10,10),doneHook='stressFree()',label='triax')
 ,NewtonIntegrator(gravity=(0,0,0),damping=0.3,label='newton')
 ]

### To achieve stress-free condition using stress control to unload the packing
def stressFree():
    print("start stress-free unloading")
    O.cell.trsf=Matrix3.Identity
    triax.stressMask=7
    triax.goal=(0,0,0)
    triax.maxStrainRate=(1,1,1)
    triax.doneHook='triaxDone()'
    triax.maxUnbalanced=10

def triaxDone():
 triax.dead=True
 O.pause()

O.run(10000000,1)

### Check total contact stress within the system
print ('Total stress (contacts) = ',getStress(O.cell.hSize[1,1]*O.cell.hSize[0,0]*O.cell.hSize[2,2]))

### Reduce the radius of all particles to further make sure there is no overlap between particles
for b in range(len(O.bodies)):
    O.bodies[b].shape.radius = O.bodies[b].shape.radius * 0.99

O.step()
print ('Total stress (contacts) after size reduction = ',getStress(O.cell.hSize[1,1]*O.cell.hSize[0,0]*O.cell.hSize[2,2]))

export.text(OUT+".isoCompression")

### Save the cell size
a = O.cell.size
np.save('cellSize',a)

MWE 2
#######################################################
######## Isotropic compression for imported dense packing ############
######################################################

from yade import pack,plot,export,ymport
import numpy as np

sp=pack.SpherePack()

O.periodic=True

RADIUS_dimension=0.0250
a = np.load('cellSize.npy')

length= a[0]
height= a[1]/3
width= a[2]
thickness=RADIUS_dimension

# friction angles
wallFRIC=90
spFRIC=0.5

# boundary conditions
PI=1.e5

### Material properties
O.materials.append(CohFrictMat(isCohesive=True,density=2500,young=1e8,poisson=0.5,frictionAngle=radians(wallFRIC),normalCohesion=-1,shearCohesion=-1,label='boxMat'))
O.materials.append(CohFrictMat(isCohesive=True,density=2500,young=1e8,poisson=0.5,frictionAngle=spFRIC,normalCohesion=-1,shearCohesion=-1,label='sphereMat'))

O.cell.hSize=Matrix3(length,0,0,0,3*height,0,0,0,width)

upBox = utils.box(center=(length/2,2*height+1*thickness,1.5*width),orientation=Quaternion(1,0,0,0),extents=(2*length,thickness/5.,width),fixed=1,wire=True,color=(1,0,0),material='boxMat')
lowBox = utils.box(center=(length/2,height-1*thickness,1.5*width),orientation=Quaternion(1,0,0,0),extents=(2*length,thickness/5.,width),fixed=1,wire=True,color=(1,0,0),material='boxMat')

O.bodies.append([upBox,lowBox])

### Import isotropic-compressed packing and boundary
packing = ymport.text("mypacking.isoCompression",color=(1,1,0),material='sphereMat')

### Filter only desired packing
def sphereWanted(b):
    x,y,z = b.state.pos
    return x >= 0 and x <= length and y >= height and y <= 2*height and z >= 0 and z <= width

packing = [b for b in packing if sphereWanted(b)]

O.bodies.append(packing)

effCellVol=(O.bodies[0].state.pos[1]-O.bodies[1].state.pos[1])*O.cell.hSize[0,0]*O.cell.hSize[2,2]
volRatio=(O.cell.hSize[0,0]*O.cell.hSize[1,1]*O.cell.hSize[2,2])/effCellVol

O.engines=[
 ForceResetter()
 ,InsertionSortCollider([Bo1_Box_Aabb(),Bo1_Sphere_Aabb()],verletDist=-0.1,allowBiggerThanPeriod=True)
 ,InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom6D()],
  [Ip2_CohFrictMat_CohFrictMat_CohFrictPhys()],
  [Law2_ScGeom6D_CohFrictPhys_CohesionMoment()]
 )
 ,GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=1,timestepSafetyCoefficient=0.8,defaultDt=-1)
 ,PeriTriaxController(dynCell=True,mass=10,maxUnbalanced=10,relStressTol=1e-4,stressMask=7,goal=(-PI/volRatio,-PI/volRatio,-PI/volRatio),globUpdate=1,maxStrainRate=(1,1,1),doneHook='triaxDone()',label='triax')
 ,NewtonIntegrator(gravity=(0,0,0),damping=0.3,label='newton')
 ]

def triaxDone():
 triax.dead=True
 O.pause()

O.step()
print ('Total stress (contacts) before isotropic compression = ',getStress((O.bodies[0].state.pos[1]-O.bodies[1].state.pos[1])*O.cell.hSize[0,0]*O.cell.hSize[2,2]))

O.run(100000000,1)

print ('Normal stress (platen) = ',O.forces.f(0)[1]/(O.cell.hSize[0,0]*O.cell.hSize[2,2]))
print ('Total stress (contacts) after isotropic compression = ',getStress((O.bodies[0].state.pos[1]-O.bodies[1].state.pos[1])*O.cell.hSize[0,0]*O.cell.hSize[2,2]))

MWE3
##########################################################
######## Isotropic compression for non-imported loose packing ############
#########################################################

from yade import pack,plot,export,ymport
import numpy as np

sp=pack.SpherePack()

O.periodic=True

RADIUS_dimension=0.0250

a = np.load('cellSize.npy')

length= a[0]
height= a[1]/3
width= a[2]
thickness=RADIUS_dimension

sp1=pack.SpherePack()

### Particle size distribution
psdSizes=[.0232,.0254,.0276]
psdCumm=[0,0.5,1]

# friction angles
wallFRIC=90
spFRIC=0.5

# boundary conditions
PI=1.e5

### Material properties
O.materials.append(CohFrictMat(isCohesive=True,density=2500,young=1e8,poisson=0.5,frictionAngle=radians(wallFRIC),normalCohesion=-1,shearCohesion=-1,label='boxMat'))
O.materials.append(CohFrictMat(isCohesive=True,density=2500,young=1e8,poisson=0.5,frictionAngle=spFRIC,normalCohesion=-1,shearCohesion=-1,label='sphereMat'))

O.cell.hSize=Matrix3(length,0,0,0,3*height,0,0,0,width)

upBox = utils.box(center=(length/2,2*height+1*thickness,1.5*width),orientation=Quaternion(1,0,0,0),extents=(2*length,thickness/5.,width),fixed=1,wire=True,color=(1,0,0),material='boxMat')
lowBox = utils.box(center=(length/2,height-1*thickness,1.5*width),orientation=Quaternion(1,0,0,0),extents=(2*length,thickness/5.,width),fixed=1,wire=True,color=(1,0,0),material='boxMat')

O.bodies.append([upBox,lowBox])

sp1.makeCloud((0,height,width),(length,2*height,2*width), psdSizes = psdSizes, psdCumm = psdCumm, periodic=True, seed =1)

sphere_id = O.bodies.append([utils.sphere(s[0],s[1],color=(1,1,0),material='sphereMat') for s in sp1])

effCellVol=(O.bodies[0].state.pos[1]-O.bodies[1].state.pos[1])*O.cell.hSize[0,0]*O.cell.hSize[2,2]
volRatio=(O.cell.hSize[0,0]*O.cell.hSize[1,1]*O.cell.hSize[2,2])/effCellVol

O.engines=[
 ForceResetter()
 ,InsertionSortCollider([Bo1_Box_Aabb(),Bo1_Sphere_Aabb()],verletDist=-0.1,allowBiggerThanPeriod=True)
 ,InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom6D()],
  [Ip2_CohFrictMat_CohFrictMat_CohFrictPhys()],
  [Law2_ScGeom6D_CohFrictPhys_CohesionMoment()]
 )
 ,GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=1,timestepSafetyCoefficient=0.8,defaultDt=-1)
 ,PeriTriaxController(dynCell=True,mass=10,maxUnbalanced=10,relStressTol=1e-4,stressMask=7,goal=(-PI/volRatio,-PI/volRatio,-PI/volRatio),globUpdate=1,maxStrainRate=(1,1,1),doneHook='triaxDone()',label='triax')
 ,NewtonIntegrator(gravity=(0,0,0),damping=0.3,label='newton')
 ]

def triaxDone():
 triax.dead=True
 O.pause()

print ('Total stress (contacts) before isotropic compression = ',getStress((O.bodies[0].state.pos[1]-O.bodies[1].state.pos[1])*O.cell.hSize[0,0]*O.cell.hSize[2,2]))

O.run(100000000,1)

print ('Normal stress (platen) = ',O.forces.f(0)[1]/(O.cell.hSize[0,0]*O.cell.hSize[2,2]))
print ('Total stress (contacts) after isotropic compression = ',getStress((O.bodies[0].state.pos[1]-O.bodies[1].state.pos[1])*O.cell.hSize[0,0]*O.cell.hSize[2,2]))

#############################################################

Best regards,
Chien-Cheng