Simulation abnormal stop “terminate called after throwing an instance of 'std::bad_alloc'” with Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM.recordMoments

Hi all,
I am trying to continue my simulation on a single particle breakage, and I used a script based on the brazilian.py[1] but changed to JCFmat materials. When I start the simulation, the simulation will stop at a certain step ( I check it with O.iter), just like falling into an endless loop. The specific specific number of steps is uncertain, sometimes occur at early stage and sometimes at middle stage, so it really bother me. Then the simulation will crush with the information:
"terminate called after throwing an instance of 'std::bad_alloc'
what(): std::bad_alloc
Aborted (core dumped)"
Maybe that is because memory is exhausted, but I can not figure out why it happened. Therefore, I want to some help.

[1]\trunk\examples\concrete

The scripts as following :
from __future__ import print_function
import math
# default parameters or from table
# generate paremters,
frictionAngle = radians(30)
# Sphere model setting
intRadius = 1.3
dtSafety = .8
strainRate = 1
specimenLength = .15
specimenRadius = .05
sphereRadius = 3.5e-3
# paramters for wall
# paramters for clump
xsize = sphereRadius * 20;
ysinze = sphereRadius * 20;
zsize = specimenRadius * 5;

# material
plate_material = JCFpmMat(young=5*200e9, density=8000, frictionAngle=frictionAngle, poisson=0.3, label='plate_mat')
sphere_material = JCFpmMat(young=4e9, cohesion=.4e8 * 0.05, density=1400, frictionAngle=frictionAngle,
                           tensileStrength=1e8 * 0.05, poisson=0.25, label='sample_mat')
O.materials.append(plate_material)
O.materials.append(sphere_material)

# spheres
pred = pack.inSphere(center=(0, 0, 0), radius=specimenRadius)
sp = pack.randomDensePack(pred, radius=sphereRadius, spheresInCell=500, memoizeDb='/tmp/packing-brazilian.db',
                          returnSpherePack=True)
sp.toSimulation(material="sample_mat")

# loading part
zMin, zMax = [pt[2] for pt in aabbExtrema()]
wallIDs = O.bodies.append([wall((0, 0, z), 2, material="plate_mat") for z in (zMin, zMax)])
walls = wallMin, wallMax = [O.bodies[i] for i in wallIDs]
v = strainRate * 2 * specimenRadius
wallMin.state.vel = (0, 0, +v)
wallMax.state.vel = (0, 0, -v)

# engines
O.engines = [
    ForceResetter(),
    InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=intRadius, label='is2aabb'), Bo1_Wall_Aabb()]),
    InteractionLoop(
        [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=intRadius, label='ss2sc'), Ig2_Wall_Sphere_ScGeom()],
        [Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(cohesiveTresholdIteration=1, label='jcf')],
        [Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(cracksFileExist=True, smoothJoint=False,
                                                           label='interactionLaw', neverErase=True, recordCracks=True,
                                                           recordMoments=True)]
    ),

    GlobalStiffnessTimeStepper(),
    # VTKRecorder(iterPeriod=iterper,initRun=True,recorders=['jcfpm','cracks','facets','moments'],label='vtk'),
    NewtonIntegrator(damping=.8, label="integrator"),
    PyRunner(iterPeriod=100, command='addPlotData()', initRun=True),
]

# plot stuff
def addPlotData():
    # forces of walls. f1 is "down", f2 is "up" (f1 needs to be negated for evlauation)
    f1, f2 = [O.forces.f(i)[2] for i in wallIDs]
    f1 *= -1
    # average force
    f = .5 * (f1 + f2)
    # displacement (2 times each wall)
    wall = O.bodies[wallIDs[0]]
    dspl = 2 * wall.state.displ()[2]
    # stress (according to standard brazilian test evaluation formula)
    stress = f / (pi * specimenRadius * specimenLength)
    # store values
    yade.plot.addData(
        t=O.time,
        i=O.iter,
        dspl=dspl,
        f1=f1,
        f2=f2,
        f=f,
        stress=stress,
    )

# plot dspl on x axis, stress on y1 axis and f,f1,f2 in y2 axis
plot.plots = {'dspl': ('stress', None, 'f1', 'f2', 'f')}

O.dt = 0.
O.step() # to create initial contacts
# now reset the interaction radius and go ahead
ss2sc.interactionDetectionFactor = 1.
is2aabb.aabbEnlargeFactor = 1.

# time step
O.dt = dtSafety * PWaveTimeStep()

# run simulation
plot.plot()
# O.run()

###
BTW, I am using this version yade-2022-03-18.git
Thanks in advance for your respopnse!