yade-users team mailing list archive

[Yuxuan Wen <question695935@xxxxxxxxxxxxxxxxxxxxx>]

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

    Status: Answered => Open

Yuxuan Wen is still having a problem:
Hello Robert and Jan,

Thank you so much for the reply!!! Now I have make a MWE version, I
delete all irrelevant and post-process commands and run the code to make
sure there is no syntax error. I wonder if it is possible that you can
help me to take a look? Thank you so much!

I found that the code can be run in Yade 2020.01a without an error, but
in Yadedaily the "Floating point exception (core dumped)" error shows
up.

The code is shown as follows, there are two files. I run the 1st file
and when it completed I run the 2nd file, since I use the O.save() at
the end of the 1st file and then use the O.load() at the beginning of
the 2nd file. Both the files are in the minimum version. The 1st file
will cost you about 1min (virtual time to 1.67s). And then when you run
the 2nd file by Yadedaily, the "Floating point exception (core dumped)"
error shows up only after a few seconds.

#-------------------------------------------------------------------------------- first file, 83 lines ---------------------------------------------------------------------------
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))

## 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) 
O.bodies.append([sphere(center,rad,material='spheres') for center,rad in sp])

######################### 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), 
	GlobalStiffnessTimeStepper(),
	NewtonIntegrator(damping=0.2), 
	PyRunner(command='finished()',iterPeriod=200),
]

######################### define functions #########################
def servo():
	O.bodies[upbox].dynamic=False
	O.bodies[lowbox].dynamic=False
	if O.cell.size[0] > (lx+1e-8):
		rate1=0.5
	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.5
	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 finished():
	if O.time>=2e4*sqrt(4.0/3.0*3.1415926*target_r*target_r*target_r*density/(kn/2)):
		O.save('consolidation.xml')
		print ('Consolidation Finished')
		O.pause()




#---------------------------------------------------------------------------------- 2nd file, 67 lines-----------------------------------------------------------------------------------
 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
ly=O.cell.size[1]

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

## record all the sphere's id
spherelist=[]
for i in O.bodies: # bodies 0 and 1 are already deleted
  spherelist.append(i.id)

## 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]

######################### 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), 
	GlobalStiffnessTimeStepper(),
	NewtonIntegrator(damping=0.2), 
	PyRunner(command='finished()',iterPeriod=1000),
]

######################### define functions #########################
def servo():
	O.bodies[upclump].dynamic=False
	O.bodies[lowclump].dynamic=False
	O.bodies[upclump].state.vel=(0,velocity,0) 	

def finished():
	if (O.time-t0)>(2*ly/velocity):
		print ('Shearing Finished')
		O.pause()

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