yade-users team mailing list archive

[Ketian Sun <question706654@xxxxxxxxxxxxxxxxxxxxx>]

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.