yade-users team mailing list archive

Thread
Date
Re: [Question #695935]: Floating point exception (core dumped)

To: yade-users@xxxxxxxxxxxxxxxxxxx
From: Yuxuan Wen <question695935@xxxxxxxxxxxxxxxxxxxxx>
Date: Tue, 09 Mar 2021 05:41:03 -0000
Reply-to: question695935@xxxxxxxxxxxxxxxxxxxxx
Sender: bounces@xxxxxxxxxxxxx
Question #695935 on Yade changed:
https://answers.launchpad.net/yade/+question/695935

Yuxuan Wen gave more information on the question:
I run this file first:
#-----------------------------------------------------------------------------------------------------
## length (m), time (s), mass (1kg), force (N), pressure (Pa)
from yade import pack, plot, qt, export, os

##
lx=0.2
ly=0.2
lz=0.14
target_r=0.005 # target r=0.005m=0.5cm, d=0.01m=1cm, l=20d
target_d=2*target_r
target_phi=0.550
target_e=(1-target_phi)/target_phi
num_spheres=int(lx*ly*lz*target_phi/(4.0/3.0*3.1415926*target_r*target_r*target_r))
flag=1 # mark time when changing friction

## space for generating particle cloud
thickness=target_d # top and bottom walls' thickness
wallsurfh=target_r # distance between wall surface to the cell
mn=Vector3(target_d,target_d,wallsurfh+thickness+target_d)
mx=Vector3(2*lx-target_d,2*ly-target_d,wallsurfh+thickness+2*lz-target_d)

## material properties
density=2650.0 # kg/m3
kn=4e5 # kn/2 is the stiffness of the contact
ks=kn*2.0/7.0
gamma=50.0
cn=4.0/3.0*3.1415926*(0.005*0.005*0.005)*density/4*gamma*2 #0.03469
frict=0 
frictnew=26.57 # after stable, change the particle from frictionless to frictional

## insert periodic cell
O.periodic=True
O.cell.hSize=Matrix3( 2*lx, 0, 0,
		      0, 2*ly, 0,
		      0, 0, wallsurfh+thickness+2*lz+thickness+wallsurfh)

## insert the walls for consolidation 
O.materials.append(ViscElMat(kn=kn,ks=ks,frictionAngle=radians(frict),cn=cn,cs=0,density=density,label='Box')) 
lowbox=O.bodies.append(utils.box(center=(lx,ly,wallsurfh+0.5*thickness), extents=(5,5,0.5*thickness), fixed=False, material='Box',wire=False)) 
upbox=O.bodies.append(utils.box(center=(lx,ly,wallsurfh+thickness+2*lz+0.5*thickness), extents=(5,5,0.5*thickness), fixed=False, material='Box',wire=False)) 

## use a SpherePack object to generate a random loose particles packing
O.materials.append(ViscElMat(kn=kn,ks=ks,frictionAngle=radians(frict),cn=cn,cs=0,density=density,label='spheres'))
sp=pack.SpherePack()
sp.makeCloud(mn,mx,rMean=target_r,num=num_spheres,periodic=True,seed=1) #"seed" make the "random" generation always the same
O.bodies.append([sphere(center,rad,material='spheres') for center,rad in sp])

msphere=0.0
sphereid=[]
for i in O.bodies:
	if i.id>=2:
		msphere=msphere+i.state.mass
		sphereid.append(i.id)


######################### define engines #########################
O.engines=[
	ForceResetter(),
	InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()],allowBiggerThanPeriod=True), 
	InteractionLoop(
		[Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()], 
		[Ip2_ViscElMat_ViscElMat_ViscElPhys()],
		[Law2_ScGeom_ViscElPhys_Basic()]
	),
	PyRunner(command='servo()',iterPeriod=1), # add force must before NewtonIntegrator(), otherwise it will be ForceResetter()
	GlobalStiffnessTimeStepper(),
	NewtonIntegrator(damping=0.2), # gravity=(0,0,0), numerical damping, dissipate Ekine because of quasi-static loading
	PyRunner(command='addPlotData()',iterPeriod=200), 
	PyRunner(command='finished()',iterPeriod=200),
	#PyRunner(command='fabricdata()',iterPeriod=500),	
	#VTKRecorder(fileName='vtkfile/consolidation-',recorders=['spheres','intr','velocity','stress','force'],iterPeriod=2000),
	#qt.SnapshotEngine(iterPeriod=2000,fileBase='video/consolidation-',label='snapshooter')
]

#plot.live=True
plot.plots={'t':('sxx','syy','szz'),'t ':('exx','eyy','ezz'),'t  ':('void'),'t   ':('Eunbal')}
plot.plot()
#O.run()

######################### define functions #########################
def servo():
	O.bodies[upbox].dynamic=False
	O.bodies[lowbox].dynamic=False
	if O.cell.size[0] > (lx+1e-8):
		rate1=0.1
	else:
		rate1=0
	O.cell.velGrad=Matrix3(-rate1,0,0, 0,-rate1,0, 0,0,0)

	if O.bodies[upbox].state.pos[2]-O.bodies[lowbox].state.pos[2]-thickness > (lz+1e-8):
		rate2=0.1
	else:
		rate2=0		
	vz=rate2*(O.bodies[upbox].state.pos[2]-O.bodies[lowbox].state.pos[2]-thickness)/2	
	O.bodies[lowbox].state.vel=(0,0,vz)
	O.bodies[upbox].state.vel=(0,0,-vz)

def addPlotData():
	V=O.cell.size[0]*O.cell.size[1]*(O.bodies[upbox].state.pos[2]-O.bodies[lowbox].state.pos[2]-thickness)
	Vs=msphere/density	
	n0=0
	n1=0
	for i in sphereid:
		if len(O.bodies[i].intrs())==0:
			n0=n0+1
		if len(O.bodies[i].intrs())==1:
			n1=n1+1

	plot.addData( 
		nsphere=num_spheres,
		N0=n0,
		N1=n1,
		step=O.iter,                
		t=O.time,	
		phi=Vs/V,
		void=(V-Vs)/Vs,	
		sxx=utils.getStress(volume=V)[0][0], sxy=utils.getStress(volume=V)[0][1], sxz=utils.getStress(volume=V)[0][2],
		syx=utils.getStress(volume=V)[1][0], syy=utils.getStress(volume=V)[1][1], syz=utils.getStress(volume=V)[1][2],
		szx=utils.getStress(volume=V)[2][0], szy=utils.getStress(volume=V)[2][1], szz=utils.getStress(volume=V)[2][2],
		exx=O.cell.size[0],
		eyy=O.cell.size[1],
		ezz=O.bodies[upbox].state.pos[2]-O.bodies[lowbox].state.pos[2]-thickness, 
		s_up=O.forces.f(upbox)[2]/(O.cell.size[0]*O.cell.size[1]),		
		s_low=O.forces.f(lowbox)[2]/(O.cell.size[0]*O.cell.size[1]),
		Z=avgNumInteractions(), # interactions between clump particle-walls are also considered in YADE 2018.02b
		Zm=avgNumInteractions(skipFree=True),
		Ekine=utils.kineticEnergy(),
		Etot=utils.O.energy.total,
		Eunbal=utils.unbalancedForce(),**O.energy
		# **O.energy converts dictionary-like O.energy to plot.addData arguments
		)
	plot.saveDataTxt(os.getcwd()+'/normalfiles/consolidation',vars=('t','step','sxx','syy','szz','exx','eyy','ezz','phi','nsphere','Eunbal','N0','N1')) # overall state variables under different step are saved in one file, no loop

def fabricdata():
	filename=os.getcwd()+'/fabricdata/'+'consolidation'+str(O.iter)
	f = open(filename,'w')
	f.write('id1 id2 x_cp y_cp z_cp n_x n_y n_z Fn_x Fn_y Fn_z Fs_x Fs_y Fs_z ovp\n')
	for i in O.interactions:
		if not i.isReal: continue
		point = i.geom.contactPoint
		norm = i.geom.normal
		ovp = i.geom.penetrationDepth
		Fn = i.phys.normalForce
		Fs = i.phys.shearForce
		f.write('%-8i %-8i %-12g %-12g %-12g %-12g %-12g %-12g %-12g %-12g %-12g %-12g %-12g %-12g %-12g\n'%(i.id1,i.id2,point[0],point[1],point[2],norm[0],norm[1],norm[2],Fn[0],Fn[1],Fn[2],Fs[0],Fs[1],Fs[2],ovp))
	f.close


def finished():
	V=O.cell.size[0]*O.cell.size[1]*(O.bodies[upbox].state.pos[2]-O.bodies[lowbox].state.pos[2]-thickness)
	Vs=msphere/density
	phi=Vs/V
	void=(V-Vs)/Vs	
	Eunbal=utils.unbalancedForce()
	global flag, step1
	if (Eunbal<=1e-6 or O.time>=1e4*sqrt(4.0/3.0*3.1415926*target_r*target_r*target_r*density/(kn/2))) and flag==1: 
		flag=2
		step1=O.iter
		setContactFriction(radians(frictnew)) # frictionAngle of boxes are still 0 as they are not dynamic body 
		print ('friction added at '+ str(O.iter))
		print ('step='+str(O.iter))
		print ('phi='+str(phi), 'error='+str(fabs((target_phi-phi)/target_phi)))
		print ('void ratio='+str(void), 'error='+str(fabs((target_e-void)/target_e)))
		print ('radius='+str(O.bodies[2].shape.radius))	
		print ('Unbalanced force = '+str(Eunbal))	
	elif (Eunbal<=1e-6 or O.time>=1e4*sqrt(4.0/3.0*3.1415926*target_r*target_r*target_r*density/(kn/2))) and O.iter-step1>=10000: 
		fabricdata()
		O.save('consolidation.xml')
		print ('Consolidation Finished')
		print ('step='+str(O.iter))
		print ('phi='+str(phi), 'error='+str(fabs((target_phi-phi)/target_phi)))
		print ('void ratio='+str(void), 'error='+str(fabs((target_e-void)/target_e)))
		print ('radius='+str(O.bodies[2].shape.radius))
		print ('Unbalanced force = '+str(Eunbal))	
		O.pause()

-- 
You received this question notification because your team yade-users is
an answer contact for Yade.