Position of particles become NaN after using Break Function

Hi everyone, I’m trying to use Break Function in a three-axis simulation, but I always get an error after doing Clump replacement: numpy.linalg.LinAlgError: Array must not contain infs or NaNs. I did some debugging and found that the error comes from some balls and Interactions having NaN values for their positions and contact points. I don’t know why they become NaN, and I wonder if anyone can help me find the reason. Here is my sample code:

from yade import pack,plot
from bf import stressTensor, checkFailure, replaceSphere, evalClump
import time

readParamsFromTable(
 key='4e5'
)
from yade.params.table import *

starttime=time.time()

unknownOk=True

#Parameters

young=8e8
poisson=0.05
density=2810
num_spheres=1000
rMean=0.18
msgoal=-10e3
msdegree=30#Friction angle during sample preparation
finalFricDegree=30#Friction angle during shearing
rate=-0.002#Shear rate (% / s)
mskr=float(key)#Rotational stiffness during sample preparation
finalkr=float(key)#Rotational stiffness during shearing
ktw=0

mn,mx=Vector3(0,0,0),Vector3(5,10,5)
stabilityThreshold=0.01

#Particle breaking control parameters
radius_ratio = 3
tension_strength=4e5
compressive_strength=4e6
wei_V0= (pi * (0.45)**3) * 4 / 3
wei_P=0.63
wei_m=2.2
discretization = 20
ifweibull=True
relative_gap = 0
grow_radius = 1
max_scale = 5
initial_packing_scale = 1.5
count_broken_particles=0
#min_radius_to_break=0.3*1e-3

#Create a material (density unit:Kg/ m ^ 3)
O.materials.append(FrictMat(young=young,poisson=poisson,frictionAngle=radians(msdegree),density=density,label='spheres'))

sp=SpherePack()

#Making Ball
sp.makeCloud(
 rMean=rMean,
 minCorner=mn,
 maxCorner=mx,
 rRelFuzz=.33,
 num=num_spheres,
 #porosity=targetPorosity,
 )
sp.toSimulation()

## Create boundaries around the Bodies
O.materials.append(FrictMat(young=young,poisson=poisson,frictionAngle=0,density=0,label='frictionless'))
walls=aabbWalls(thickness=1e-10,material='frictionless')
wallIds=O.bodies.append(walls)

triax=TriaxialStressController(
 wall_bottom_id=wallIds[2],
 wall_top_id=wallIds[3],
 wall_left_id=wallIds[0],
 wall_right_id=wallIds[1],
 wall_back_id=wallIds[4],
 wall_front_id=wallIds[5],
 internalCompaction=True,
 stressMask=7,
 goal1=msgoal,
 goal2=msgoal,
 goal3=msgoal,
 maxMultiplier=1. + 1e-4, # spheres growing factor (fast growth)
 finalMaxMultiplier=1. + 1e-5, # spheres growing factor (slow growth)
 max_vel=1,#（m/s）
 label="triax",
 #thickness=0,
)

cbcontroller=PyRunner(iterPeriod=10000,command='clumpUpdate(radius_ratio, tension_strength, compressive_strength, wei_V0, wei_P,wei_m)',label='clumpbreakagecontroller')

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_MindlinPhys(
   krot=mskr,
   ktwist=ktw,
   )
  ],
  [Law2_ScGeom_MindlinPhys_Mindlin(
   includeMoment=True,
  )]
 ),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
 triax,
 PyRunner(iterPeriod=200,command='history()',label='triaxrecorder'),
 cbcontroller,
 NewtonIntegrator(damping=.3)
 ]
O.trackEnergy = True

## a function saving variables
def history():
  plot.addData(
   e11=-triax.strain[0]*100, e22=-triax.strain[1]*100, e33=-triax.strain[2]*100,
   ev=100*(triax.strain[0]+triax.strain[1]+triax.strain[2]),
   s11=-triax.stress(0)[0]/1000,
   s22=-triax.stress(2)[1]/1000,
   s33=-triax.stress(4)[2]/1000,
   diffstress=(-triax.stress(2)[1]-(-triax.stress(0)[0]-triax.stress(4)[2])/2)/1000,
   Etot=O.energy.total(),
   i=O.iter,
   **O.energy
   )
plot.plots={
 'i':('e11','e22','e33',),
 'i ':('s11','s22','s33'),
 ' i':(O.energy.keys,None,'Etot'),
 'e22':('s11','s22','s33',None,'ev'),
 }

def updateFile(file):
 """
 mod record txt for OrigionLab
 """
 import re
 f=open(file, "r", encoding="utf-8")
 l=f.readlines()
 l[0]=re.sub(r'[#\t]\t','\t',l[0])
 del l[1]
 f=open(file,'w',encoding="utf-8")
 f.writelines(l)
 f.close()

def clumpUpdate(radius_ratio, tension_strength, compressive_strength, wei_V0, wei_P,wei_m):
 mn_box=Vector3(O.bodies[wallIds[0]].state.pos[0],O.bodies[wallIds[2]].state.pos[1],O.bodies[wallIds[4]].state.pos[2])
 mx_box=Vector3(O.bodies[wallIds[1]].state.pos[0],O.bodies[wallIds[3]].state.pos[1],O.bodies[wallIds[5]].state.pos[2])
 outer_predicate=pack.inAlignedBox(mn_box,mx_box)

 global count_broken_particles
 for i in range(
   len(O.bodies)
 ):
  b = O.bodies[i]
  if not b == None:
   if isinstance(b.shape, Clump):
    clump_broken = evalClump(
      b.id,
      radius_ratio,
      tension_strength,
      compressive_strength,
      outer_predicate=outer_predicate,
      max_scale=max_scale,
      grow_radius=grow_radius,
      relative_gap=relative_gap,
      discretization=discretization,
      weibull=True,
      wei_V0=wei_V0,
      wei_P=wei_P,
      wei_m=wei_m
    )
    # refresh time step
    if clump_broken:
     count_broken_particles += 1
 print('Broken particles:{}'.format(count_broken_particles))

def mainlogic():
 #stage 1
 cbcontroller.dead=True
 while 1:
  O.run(1000, True)
  unb=unbalancedForce()

  print('Stage 1 :unbalanced force:{:.3},mean Stress:{:.3},porosity:{:.4},speed:{:d},current step:{:d}'.format(unb,triax.meanStress,triax.porosity,int(O.speed),O.iter))
  if unb<stabilityThreshold and abs(msgoal-triax.meanStress)/(-msgoal)<0.001:
   print("### Stage 1 Complete! Balls Created ###")
   break
 midtime=time.time()
 print('time spended:'+str(midtime-starttime)+'s')

 relRadList2 = [.5, 1, .5]
 relPosList2 = [[1, 0, 0], [0, 0, 0], [-1, 0, 0]]
 templates = []
 templates.append(clumpTemplate(relRadii=relRadList2, relPositions=relPosList2))
 #replace by 10%
 O.bodies.replaceByClumps(templates, [0.1],discretization=20)

 triax.goal2=msgoal
 triax.goal1=msgoal
 triax.goal3=msgoal

 #stage 2
 while 1:
  O.run(1000, True)
  unb=unbalancedForce()
  print('Stage 2 :unbalanced force:{:.3},mean Stress:{:.3},porosity:{:.4},speed:{:d},current step:{:d}'.format(unb,triax.meanStress,triax.porosity,int(O.speed),int(O.iter)))
  if unb<stabilityThreshold and abs(msgoal-triax.meanStress)/(-msgoal)<0.001:
   print("### Stage 2 Complete! Clumps Created ###")
   break
 print('time spended:'+str(time.time()-midtime)+'s')
 midtime=time.time()

 cbcontroller.dead=False

 #Isotropic loading

 compressgoal=-400e3
 triax.internalCompaction=False
 triax.goal2=compressgoal
 triax.goal1=compressgoal
 triax.goal3=compressgoal
 compFricDegree=0
 setContactFriction(radians(compFricDegree))
 while 1:
  O.run(1000, True)
  #the global unbalanced force on dynamic bodies, thus excluding boundaries, which are not at equilibrium
  unb=unbalancedForce()
  print('Isotropic loading:unbalanced force:{:.3},mean Stress:{:.3},porosity:{:.4},speed:{:d},current step:{:d}'.format(unb,triax.meanStress,triax.porosity,int(O.speed),int(O.iter)))
  if unb<stabilityThreshold and abs(compressgoal-triax.meanStress)/(-compressgoal)<0.001:
   break
 print("### Isotropic loading Complete ###")

 triax.internalCompaction=False

 setContactFriction(radians(finalFricDegree))
 plot.plot()
 #set stress control on x and z, we will impose strain rate on y
 triax.stressMask = 5
 #now goal2 is the target strain rate
 triax.goal2=rate
 # we define the lateral stresses during the test, here the same 10kPa as for the initial confinement.
 triax.goal1=compressgoal
 triax.goal3=compressgoal
 triax.height0=triax.height
 triax.width0=triax.width
 triax.depth0=triax.depth

 for i in O.interactions:
  i.phys.kr=finalkr
 while 1:
  O.run(1000, True)
  print('Partial stress:',-triax.stress(2)[1]+triax.stress(0)[0],'Body strain',-triax.strain[0]-triax.strain[1]-triax.strain[2],'Axial strain',-triax.strain[1],'Speed:',O.speed)
  if -triax.strain[1]>0.22:
   break
 plot.saveDataTxt('data'+key+'.txt')
 updateFile('data'+key+'.txt')

mainlogic()

I have searched online for similar problems, but I couldn’t find any solution. Any help or suggestion would be greatly appreciated. Thank you in advance.