yade-users team mailing list archive

[Yuxuan Wen <question695935@xxxxxxxxxxxxxxxxxxxxx>]

Question #695935 on Yade changed:
https://answers.launchpad.net/yade/+question/695935

Yuxuan Wen gave more information on the question:
And then I run the second file. The "
Floating point exception (core dumped)" error will show when running this file:

#---------------------------------------------------------------------
## length (m), time (s), mass (1kg), force (N), pressure (Pa)
from yade import pack, plot, qt, export, os

O.load('consolidation.xml')
step0=O.iter
t0=O.time
rate=0.5
velocity=rate*0.1
radius=O.bodies[2].shape.radius
diameter=2*radius

## delete 2 up and low boxes
z1=O.bodies[0].state.pos[2] 
z2=O.bodies[1].state.pos[2]
thickness=diameter 
upzsurf=max(z1,z2)-0.5*thickness
lowzsurf=min(z1,z2)+0.5*thickness
lx=O.cell.size[0]
ly=O.cell.size[1]
lz=upzsurf-lowzsurf
V=lx*ly*lz
Vs0=O.cell.size[0]*O.cell.size[1]*O.cell.size[2]*(1-porosity())

O.bodies.erase(0,True) 
O.bodies.erase(1,True)

## record all the sphere's id
spherelist=[]
msphere=0.0
for i in O.bodies: # bodies 0 and 1 are already deleted
  spherelist.append(i.id)
  msphere=msphere+i.state.mass
Vs=msphere/2650.0

## clump the walls 
uplist=[]
lowlist=[]
clumplist=[]
for i in O.bodies:
  if i.state.pos[2] >= (upzsurf-1.5*diameter):
    uplist.append(i.id)
  if i.state.pos[2] <= (lowzsurf+1.5*diameter):
    lowlist.append(i.id)
clumplist=uplist+lowlist
upclump=O.bodies.clump(uplist)
lowclump=O.bodies.clump(lowlist)

upposy0=O.bodies[upclump].state.pos[1]
lowposy0=O.bodies[lowclump].state.pos[1]

## record the data of sphere body that is not in the clump walls
shearspherelist=[]
checklist=[]
for i in spherelist:
  if i not in uplist and i not in lowlist:
    shearspherelist.append(i)	
    checklist.append(i)	
checklist.append(upclump)
checklist.append(lowclump)


#O.bodies[upclump].dynamic=False
#O.bodies[lowclump].dynamic=False
#O.bodies[upclump].state.vel=(0,velocity,0)
######################### define engines #########################
O.engines=[
	ForceResetter(),
	InsertionSortCollider([Bo1_Sphere_Aabb()]), 
	InteractionLoop(
		[Ig2_Sphere_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(gravity=(0,0,0),damping=0.2), # numerical damping, dissipate Ekine because of quasi-static loading
	PyRunner(command='if O.iter-step0<=10000:addPlotData()',iterPeriod=100), 
	PyRunner(command='if O.iter-step0>10000:addPlotData()',iterPeriod=1000), 
	PyRunner(command='if O.iter-step0<=100000:veldata()',iterPeriod=1000), 
	PyRunner(command='if O.iter-step0<=10000:fabricdata()',iterPeriod=100),
	PyRunner(command='if O.iter-step0>10000:fabricdata()',iterPeriod=1000),
	PyRunner(command='finished()',iterPeriod=1000),
]
plot.plots={'t':('sxx','syy','szz'),'t ':('upfx','upfy','upfz'),'t  ':('upvx','upvy','upvz'),'t   ':('Eunbal')}
plot.plot()

######################### define functions #########################
def servo():
	O.bodies[upclump].dynamic=False
	O.bodies[lowclump].dynamic=False
	O.bodies[upclump].state.vel=(0,velocity,0) # can't blockDOF and add velocity at the same time when dynamic is turned to False

def addPlotData():
	upforcex=0
	upforcey=0
	upforcez=0
	lowforcex=0
	lowforcey=0
	lowforcez=0
	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
		if (i.id1 in uplist and i.id2 not in uplist): # force direction is from id1 to id2  
			upforcex=upforcex-Fn[0]-Fs[0]
			upforcey=upforcey-Fn[1]-Fs[1]
			upforcez=upforcez-Fn[2]-Fs[2]
		if (i.id1 not in uplist and i.id2 in uplist): 
			upforcex=upforcex+Fn[0]+Fs[0]	
			upforcey=upforcey+Fn[1]+Fs[1]
			upforcez=upforcez+Fn[2]+Fs[2]
		if (i.id1 in lowlist and i.id2 not in lowlist):
			lowforcex=lowforcex-Fn[0]-Fs[0]	
			lowforcey=lowforcey-Fn[1]-Fs[1]
			lowforcez=lowforcez-Fn[2]-Fs[2]
		if (i.id1 not in lowlist and i.id2 in lowlist):
			lowforcex=lowforcex+Fn[0]+Fs[0]	
			lowforcey=lowforcey+Fn[1]+Fs[1]
			lowforcez=lowforcez+Fn[2]+Fs[2]
	n0=0
	n1=0
	for i in shearspherelist:
		if len(O.bodies[i].intrs())==0:
			n0=n0+1
		if len(O.bodies[i].intrs())==1:
			n1=n1+1

	plot.addData( 
		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.cell.size[2], 
		upfx=upforcex, upfy=upforcey, upfz=upforcez,
		lowfx=lowforcex, lowfy=lowforcey, lowfz=lowforcez,		
		upclumpfx=O.forces.f(upclump)[0], upclumpfy=O.forces.f(upclump)[1], upclumpfz=O.forces.f(upclump)[2],
		lowclumpfx=O.forces.f(lowclump)[0], lowclumpfy=O.forces.f(lowclump)[1], lowclumpfz=O.forces.f(lowclump)[2],		
		upvx=O.bodies[upclump].state.vel[0], upvy=O.bodies[upclump].state.vel[1], upvz=O.bodies[upclump].state.vel[2], 
		lowvx=O.bodies[lowclump].state.vel[0],lowvy=O.bodies[lowclump].state.vel[1],lowvz=O.bodies[lowclump].state.vel[2],
		upposx=O.bodies[upclump].state.pos[0],upposy=O.bodies[upclump].state.pos[1],upposz=O.bodies[upclump].state.pos[2],
		lowposx=O.bodies[lowclump].state.pos[0],lowposy=O.bodies[lowclump].state.pos[1],lowposz=O.bodies[lowclump].state.pos[2],
		Z=avgNumInteractions(), # interactions between clump particle-walls are also considered in YADE 2018.02b
		Zm=avgNumInteractions(skipFree=True),
		nshearsphere=len(shearspherelist), nclumpsphere=len(clumplist), nsphere=len(spherelist),
		Ekine=utils.kineticEnergy(),
		Etot=utils.O.energy.total,
		Eunbal=utils.unbalancedForce(),**O.energy # **O.energy converts dictionary-like O.energy to plot.addData arguments		
		)
	# overall state variables under different step are saved in one file, no loop
	plot.saveDataTxt(os.getcwd()+'/normalfiles/shearstate',vars=('t','step','Eunbal','phi','void','nshearsphere','nclumpsphere','nsphere')) 
	plot.saveDataTxt(os.getcwd()+'/normalfiles/shearforce',vars=('t','step','upfx','upfy','upfz','lowfx','lowfy','lowfz','upposy','upvy','lowposy','lowvy','N0','N1'))
	#plot.saveDataTxt(os.getcwd()+'/normalfiles/shearflowclump',vars=('t','step','upclumpfx','upclumpfy','upclumpfz','lowclumpfx','lowclumpfy','lowclumpfz','upvx','upvy','upvz','lowvx','lowvy','lowvz'))


def veldata():
	filename=os.getcwd()+'/velocitydata/'+'velocity'+str(O.iter)
	f = open(filename,'w')
	f.write('id vx vy vz pos_x pos_y pos_z ovp\n')
	for i in checklist: # data of all bodies/intrs at each step is saved under a unique file named str(O.iter) 
		vx=O.bodies[i].state.vel[0]
		vy=O.bodies[i].state.vel[1]
		vz=O.bodies[i].state.vel[2]
		pos_x=O.bodies[i].state.pos[0]
		pos_y=O.bodies[i].state.pos[1]
		pos_z=O.bodies[i].state.pos[2]
		f.write('%-8i %-12g %-12g %-12g %-12g %-12g %-12g\n'%(i,vx,vy,vz,pos_x,pos_y,pos_z))
	f.close

def fabricdata():
	filename1=os.getcwd()+'/fabricdata/'+'fabric'+str(O.iter)
	f1 = open(filename1,'w')
	f1.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
		if not (i.id1 in clumplist and i.id2 in clumplist): # the clumplist interactions should be constant 			
			f1.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))
	f1.close	
		
	filename2=os.getcwd()+'/N0N1data/'+'N0_id'+str(O.iter)
	f2 = open(filename2,'w')
	f2.write('N0_id\n')
	for i in shearspherelist:
		if len(O.bodies[i].intrs())==0:
			f2.write('%-8i\n'%(i))
	f2.close
	
	filename3=os.getcwd()+'/N0N1data/'+'N1_id'+str(O.iter)
	f3 = open(filename3,'w')
	f3.write('N1_id\n')
	for i in shearspherelist:
		if len(O.bodies[i].intrs())==1:
			f3.write('%-8i\n'%(i))
	f3.close

def finished():
	unbal=utils.unbalancedForce()
	upvely=O.bodies[upclump].state.vel[1]
	upposy=O.bodies[upclump].state.pos[1]
	strain=(upposy-upposy0)/lz

	if (O.time-t0)>(20*ly/velocity):
		print ('Shearing Finished')
		print ('Unbalanced force='+str(unbal))	
		print ('Shear strain='+str(strain))
		O.pause()

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