Floating point exception (core dumped)

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()

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()

Yikes, that's a lot of files. Any chance you can make an MWE [1]?

[1]https://www.yade-dem.org/wiki/Howtoask

> I was using Ubuntu 18.04 ... my Ubuntu updated to 20.04 automatically

Maybe "accidentally" instead of "automatically? ("automatic upgrade" does not sound like how Ubuntu (used to?) work)

> If I want to make the code works, should I still install the Yade 2018.02b?

definitely it is an option.
Depending on how you use the Ubuntu (I have there all job and personal data, program settings, ..., would not want to install it from scratch)

> in Yade 2018.02b ... in Yade 2020.01

you have 2 years of development, some behavior may change. Concerning makeCloud and seed, there were some bugs / bug fixing [2].
I guess a significant change is Python 2 -> Python 3, which also may break something in your scripts (e.g. int/float division).

> And then I run the second file. The "Floating point exception (core dumped)" error will show when running this file:

what is a typical running time of the scripts?
"a lot of files" is not a big problem for mere testing.
Of course, debugging and searching for problem solution is completely different topic, where a MWE required by Robert would/should help significantly.

Concerning seed in different versions, see the commit [2], the related bug, date of the commit w.r.t yade versions you use, possibly other history of this part of the file, ......

cheers
Jan

[2] https://gitlab.com/yade-dev/trunk/-/commit/cab39356c5b3bec1473b68fdf0b8579d804f03ca

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()

I do not have 20.04, so I am out of testing, sorry..

Hello Jan,

The code can be run in Yade 2020.1a, but can not be run correctly in Yadedaily. I wonder if it is possible that you could help me to take a look on Yadedaily. It will not take you much time, for the 1st file is about 1 min, and for the 2nd file is about several seconds.
I know this is a troublesome request and will make you inconvenient, I am deeply sorry!! But I really have no clue on how to solve the issue or debugging. If you could help me to take a look or just give me some clues to solve this problem, I will be much appreciated!! Otherwise I have no choice but to install the Ubuntu 18.04 and Yade 2018.02b again...

Kind Regards,
Yuxuan

On Ubuntu 18.04 using yadedaily (Yade 20210217-4980~0ed545f~bionic1), I did not get "Floating point exception"...
Jan

Thank you so much Jan!! Your comment has strengthen my confidence on Yadedaily. I am using Yadedaily in Ubuntu 20.04 and the "Floating point exception (core dumped)" shows up. So it seems that the error is due to the version of Ubuntu instead of Yadedaily. I will install the Ubuntu 18.04 rightaway, and the problem can be solved.

Thank you so much for you time and patience!

Best Regards,
Yuxuan

Hi,
I confirm segfault on ubuntu20/yadedaily (20210204-4960~ab01ee8~focal1)
The problem is apparently with InsertionSortCollider at iteration 2 in script2.
I suspect the combination of reloading/erasing/periodicity does wrong at some point.

> I will install the Ubuntu 18.04

Oh no... don't do that. Why not Yade 2020.1a if it works for you?

It would be helpful if you could reproduce in a single script where we wouldn't have to click play to see the error (and possibly taking less than minutes to run).

Regards

Bruno

Hi,
The crash was due to erasing just after reloading.
The problem is fixed in source code in [1]. The fix is not yet in yadedaily but it should appear quickly.
If you use a yade version affected by this bug you can workaround with "O.bodies.enableRedirection=False".
Regards
Bruno

[1] https://gitlab.com/yade-dev/trunk/-/merge_requests/637

Hello Dr. Chareyre,

Thank you so much for your answer and your time on the debugging! I am sorry that I didn't reproduce a single script and reply to you last week. I am devoting in producing the results of simulation and the post-processing for a PPT so that I forgot to produce this single script. I am sorry for the negligence! It is so delightful that I can get your answer and guidance on an problem that I met and couldn't solve! Thank you for your answer and it is so helpful for me to get a deeper understanding on Yade by tracking this bug and debugging process!

Best Regards,
Yuxuan