New question #706654 on Yade:
https://answers.launchpad.net/yade/+question/706654
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.
--
You received this question notification because your team yade-users is
an answer contact for Yade.
_______________________________________________
Mailing list: https://launchpad.net/~yade-users
Post to : [email protected]
Unsubscribe : https://launchpad.net/~yade-users
More help : https://help.launchpad.net/ListHelp
