Particles fall out the container

import random
import math
from yade import geom, pack, utils, ymport
from yade import export
import numpy as np

# Define cylinder parameters
center = (0, 0, 0)
radius = 0.102
height = 0.064

# create cylindrical body with radius 0.102 m and height 0.064 m
cylinder = yade.geom.facetCylinder(center=center, radius=radius, height=height, segmentsNumber=200, wallMask=6)

# define material properties
mat = PolyhedraMat()
mat.density = 2600 # kg/m^3
mat.young = 1E6 # Pa
mat.poisson = 20000 / 1E6
mat.frictionAngle = 0.6 # rad

# assign material to each body in the cylinder
for body in cylinder:
    body.bodyMat = mat

# add cylinder to simulation
O.bodies.append(cylinder)

# Generate sphere pack

O.materials.append(FrictMat(young=1E6, poisson=20000/1E6, density=10, label='agg'))

sp = pack.SpherePack()
sp.makeCloud((-0.05,-0.05,0), (0.05,0.05,0.09),rMean=0.01575,rRelFuzz=0,periodic=False,num=6)

sp2_1 = pack.SpherePack()
sp2_1.makeCloud((-0.03,-0.03,0.1), (0.03,0.03,0.15),rMean=0.01175,rRelFuzz=0,periodic=False,num=7)

sp2_2 = pack.SpherePack()
sp2_2.makeCloud((-0.03,-0.03,0.16), (0.03,0.03,0.2),rMean=0.011,rRelFuzz=0,periodic=False,num=8)

sp2_3 = pack.SpherePack()
sp2_3.makeCloud((-0.03,-0.03,0.21), (0.03,0.03,0.25),rMean=0.01025,rRelFuzz=0,periodic=False,num=8)

sp3_1 = pack.SpherePack()
sp3_1.makeCloud((-0.03,-0.03,0.26), (0.03,0.03,0.3),rMean=0.0083125,rRelFuzz=0,periodic=False,num=31)

sp3_2 = pack.SpherePack()
sp3_2.makeCloud((-0.03,-0.03,0.31), (0.03,0.03,0.35),rMean=0.007125,rRelFuzz=0,periodic=False,num=31)

sp3_3 = pack.SpherePack()
sp3_3.makeCloud((-0.03,-0.03,0.36), (0.03,0.03,0.4),rMean=0.0059375,rRelFuzz=0,periodic=False,num=32)

sp4_1 = pack.SpherePack()
sp4_1.makeCloud((-0.03,-0.03,0.41), (0.03,0.03,0.45),rMean=0.0041525,rRelFuzz=0,periodic=False,num=25)

sp4_2 = pack.SpherePack()
sp4_2.makeCloud((-0.03,-0.03,0.46), (0.03,0.03,0.5),rMean=0.003555,rRelFuzz=0,periodic=False,num=26)

sp4_3 = pack.SpherePack()
sp4_3.makeCloud((-0.03,-0.03,0.51), (0.03,0.03,0.55),rMean=0.0029575,rRelFuzz=0,periodic=False,num=26)

sp5_1 = pack.SpherePack()
sp5_1.makeCloud((-0.03,-0.03,0.56), (0.03,0.03,0.6),rMean=0.002065,rRelFuzz=0,periodic=False,num=25)

sp5_2 = pack.SpherePack()
sp5_2.makeCloud((-0.03,-0.03,0.61), (0.03,0.03,0.65),rMean=0.00177,rRelFuzz=0,periodic=False,num=26)

sp5_3 = pack.SpherePack()
sp5_3.makeCloud((-0.03,-0.03,0.66), (0.03,0.03,0.7),rMean=0.001475,rRelFuzz=0,periodic=False,num=26)

sp6_1 = pack.SpherePack()
sp6_1.makeCloud((-0.03,-0.03,0.56), (0.03,0.03,0.6),rMean=0.001035,rRelFuzz=0,periodic=False,num=155)

sp6_2 = pack.SpherePack()
sp6_2.makeCloud((-0.03,-0.03,0.61), (0.03,0.03,0.65),rMean=0.00089,rRelFuzz=0,periodic=False,num=155)

sp6_3 = pack.SpherePack()
sp6_3.makeCloud((-0.03,-0.03,0.66), (0.03,0.03,0.7),rMean=0.000745,rRelFuzz=0,periodic=False,num=156)

sp7_1 = pack.SpherePack()
sp7_1.makeCloud((-0.03,-0.03,0.71), (0.03,0.03,0.75),rMean=0.000525,rRelFuzz=0,periodic=False,num=155)

sp7_2 = pack.SpherePack()
sp7_2.makeCloud((-0.03,-0.03,0.76), (0.03,0.03,0.8),rMean=0.00045,rRelFuzz=0,periodic=False,num=155)

sp7_3 = pack.SpherePack()
sp7_3.makeCloud((-0.03,-0.03,0.81), (0.03,0.03,0.85),rMean=0.000375,rRelFuzz=0,periodic=False,num=156)

sp8_1 = pack.SpherePack()
sp8_1.makeCloud((-0.03,-0.03,0.86), (0.03,0.03,0.9),rMean=0.0002625,rRelFuzz=0,periodic=False,num=155)

sp8_2 = pack.SpherePack()
sp8_2.makeCloud((-0.03,-0.03,0.91), (0.03,0.03,0.95),rMean=0.000225,rRelFuzz=0,periodic=False,num=155)

sp8_3 = pack.SpherePack()
sp8_3.makeCloud((-0.03,-0.03,0.96), (0.03,0.03,1),rMean=0.0001875,rRelFuzz=0,periodic=False,num=156)

sp9 = pack.SpherePack()
sp9.makeCloud((-0.03,-0.03,0.01), (0.03,0.03,0.05),rMean=0.000075,rRelFuzz=0,periodic=False,num=8000)

sp10 = pack.SpherePack()
sp10.makeCloud((-0.03,-0.03,0.04), (0.03,0.03,1),rMean=0.00005,rRelFuzz=0,periodic=False,num=80000)

# add the sphere pack to the simulation
sp.toSimulation(material='agg')
sp2_1.toSimulation(material='agg')
sp2_2.toSimulation(material='agg')
sp2_3.toSimulation(material='agg')
sp3_1.toSimulation(material='agg')
sp3_2.toSimulation(material='agg')
sp3_3.toSimulation(material='agg')
sp4_1.toSimulation(material='agg')
sp4_2.toSimulation(material='agg')
sp4_3.toSimulation(material='agg')
sp5_1.toSimulation(material='agg')
sp5_2.toSimulation(material='agg')
sp5_3.toSimulation(material='agg')
sp6_1.toSimulation(material='agg')
sp6_2.toSimulation(material='agg')
sp6_3.toSimulation(material='agg')
sp7_1.toSimulation(material='agg')
sp7_2.toSimulation(material='agg')
sp7_3.toSimulation(material='agg')
sp8_1.toSimulation(material='agg')
sp8_2.toSimulation(material='agg')
sp8_3.toSimulation(material='agg')
sp9.toSimulation(material='agg')
sp10.toSimulation(material='agg')

# Define gravity engine
O.engines = [ ForceResetter(), InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()]),
   InteractionLoop(
       [Ig2_Sphere_Sphere_ScGeom6D(), Ig2_Facet_Sphere_ScGeom()],
       [Ip2_FrictMat_FrictMat_FrictPhys()],
       [Law2_ScGeom_FrictPhys_CundallStrack()]
   ),
   GlobalStiffnessTimeStepper(defaultDt=0.0001, timestepSafetyCoefficient=0.2),
   NewtonIntegrator(gravity=(0,0,-1000), damping=0.4),
]

# Define simulation duration and time step
O.dt = 0.5e-6
O.run(10, wait=True)
for body in O.bodies:
   if not isinstance(body.shape, Sphere):
       continue
   if body.shape.radius == 0.01575: #SP
       body.shape.color = (0,0,0)

   if body.shape.radius == 0.01175 :
       body.shape.color = (0,0,1)
   if body.shape.radius == 0.011 :
       body.shape.color = (0,0,1)
   if body.shape.radius == 0.01025 :
       body.shape.color = (0,0,1)

   if body.shape.radius == 0.00831250 : #SP3
       body.shape.color = (0,1,0)
   if body.shape.radius == 0.007125 :
       body.shape.color = (0,1,0)
   if body.shape.radius == 0.0059375 :
       body.shape.color = (0,1,0)

   if body.shape.radius == 0.0041525: #SP4
       body.shape.color = (1,0,0)
   if body.shape.radius == 0.003555 :
       body.shape.color = (1,0,0)
   if body.shape.radius == 0.0029575 :
       body.shape.color = (1,0,0)

   if body.shape.radius == 0.002065: #SP5
       body.shape.color = (1,1,0)
   if body.shape.radius == 0.00177 :
       body.shape.color = (1,1,0)
   if body.shape.radius == 0.001475 :
       body.shape.color = (1,1,0)

   if body.shape.radius == 0.001035: #SP6
       body.shape.color = (1,1,1)
   if body.shape.radius == 0.00089 :
       body.shape.color = (1,1,1)
   if body.shape.radius == 0.000745 :
       body.shape.color = (1,1,1)

   if body.shape.radius == 0.000525: #SP7
       body.shape.color = (1,0.5,1)
   if body.shape.radius == 0.00045 :
       body.shape.color = (1,0.5,1)
   if body.shape.radius == 0.000375 :
       body.shape.color = (1,0.5,1)

   if body.shape.radius == 0.0002625: #SP8
       body.shape.color = (1,1,0.5)
   if body.shape.radius == 0.000225 :
       body.shape.color = (1,1,0.5)
   if body.shape.radius == 0.0001875 :
       body.shape.color = (1,1,0.5)

   if body.shape.radius == 0.000075: #SP9
       body.shape.color = (0.5,1,0.5)

   if body.shape.radius == 0.00005: #SP10
       body.shape.color = (0.5,0.5,1)

def checkUnbalanced():
 if unbalancedForce() < .05:
  O.pause()
  plot.saveDataTxt('bbb.txt.bz2')
  # plot.saveGnuplot('bbb') is also possible

# Count the number of spheres in the cylinder

num_spheres_in_cylinder = 0
for body in O.bodies:
   if body.id == 0: # skip the cylinder
       continue
   sphere_center = body.state.pos
   dist_to_axis = math.sqrt(sphere_center[0]**2 + sphere_center[1]**2)
   if dist_to_axis <= radius and sphere_center[2] >= center[2] and sphere_center[2] <= center[2] + height:
       num_spheres_in_cylinder += 1

# Calculate volume of spheres in cylinder
volume_of_spheres_in_cylinder = 0
for body in O.bodies:
   sphere_center = body.state.pos
   dist_to_axis = math.sqrt(sphere_center[0]**2 + sphere_center[1]**2)
   if dist_to_axis <= radius and sphere_center[2] >= center[2] and sphere_center[2] <= center[2] + height:
       if isinstance(body.shape, Sphere): # check if body is a sphere
           sphere_volume = (4/3) * math.pi * body.shape.radius**3
           volume_of_spheres_in_cylinder += sphere_volume

# Calculate volume of cylinder
volume_of_cylinder = (math.pi * radius**2) * height

# Calculate porosity
porosity = (volume_of_cylinder - volume_of_spheres_in_cylinder) / volume_of_cylinder
porosity_percent = porosity * 100

# Print results
print("Number of spheres in cylinder:", num_spheres_in_cylinder)
print("Volume of spheres in cylinder:", volume_of_spheres_in_cylinder)
print("Volume of cylinder:", volume_of_cylinder)
print("Porosity:", porosity)
print("Porosity:", porosity_percent, "%")

When I try to simulate I don't understand why the particles fall under the container?