Model reproducibility

Hi Andy,

I haven't tried running your script yet but looks like your verlet distance might be a tad large. From the geometry file, your range of particle radii is 0.4 to 1.0. Best to try a verlet distance <= 0.2*Rmin = 0.08 rather than 0.1 as currently implemented. Your timestep increment (0.03) appears to be small enough.

Let me know how you go with a smaller verlet distance. I'm not sure if this will fix the reproducibility problem though...

Cheers,

Dion

Very strange stuff indeed. The Verlet distance IMHO shouldn't be an issue as long as it is smaller than grid_spacing-2*r_max , which it is in your case.
Btw, what version of esys-particle is this run with? There was a potentially nasty bug in the circular boundary condition code fixed in rev. 1010.
Anyway, I'll have a shot at your script when I find the time in the next couple of days.

Steffen

I am using 2.1, which was reinstalled and last updated 8 Sept, so probably update 1008.

I am trying with a different verlet distance now and changed some other things. I can post a image tomorrow when they are done.

AndyR

Steffen, Dion (and anyone else)

I changed some parameters and get exact replication of my runs. I added interactions between the bonded walls and frictional particles and some other things. I also added a time where I only apply normal stress at the start and no shear displacement for 5000 timesteps. I don't know if this or something else fixed it (missing interaction or exclusion). I am running more repeat tests now to check. If I subtract the friction measurements from two runs I get a vector of zeros.

Thanks for any time you spent on this.
AndyR

My new code is pasted below.

from esys.lsm import*
from esys.lsm.util import*
from time import*

#define a Runnable class for compression loading during the simulation
class Loading(Runnable):
 def __init__(self,sim,dy,t0,fn):
  Runnable.__init__(self)
  self.sim=sim
  self.dy=dy
  self.counter=0
  self.t0=t0
  self.fn=fn
  self.start_time=time()

 def run(self):
  if(self.counter>self.t0+5000):
   self.sim.moveWallBy("lowerWall",Vec3(self.dy,0.0,0.0))
   self.sim.moveWallBy("upperWall",Vec3(-1.0*self.dy,0.0,0.0))
  elif(self.counter<=5000): #build dead space in start of model for compression
   self.sim.moveWallBy("lowerWall",Vec3(0.0,0.0,0.0))
   self.sim.moveWallBy("upperWall",Vec3(0.0,0.0,0.0))
  else: # this part does the ramp
   dyt=(float(self.counter-5000)/float(self.t0))*self.dy
   self.sim.moveWallBy("lowerWall",Vec3(dyt,0.0,0.0))
   self.sim.moveWallBy("upperWall",Vec3(-1.0*dyt,0.0,0.0))
  self.sim.applyForceToWall("lowerWallInteraction",Vec3(0.0,self.fn,0.0))
  self.sim.applyForceToWall("upperWallInteraction",Vec3(0.0,-1.0*self.fn,0.0))

  self.counter=self.counter+1
  if((self.counter%100)==0):
   runtime=time()-self.start_time
   npart=self.sim.getNumParticles()
   print "timestep: ",self.counter,self.counter," at ", runtime, " Total Particles: ", npart

def runSimulation():

 cpux=3
 cpuy=2
 cpuz=1
 geometry_name="/home/rathbuna/data/geometry/scripts/teeth_markers.geo"
 t0=5000 #the time we want to be at for full speed
 nt=1000000
 dt=0.015
 dt_snap=10000
 dt_force=10000
 dt_save=50000
 xmax=60.0
 ymax=58.0
 zmax=18.0
 dy=0.000015 #change in y each time step
 fn=0.0005*xmax*zmax #area times stress.

 shear_mix=LsmMpi(cpux*cpuy*cpuz,[cpux,cpuy,cpuz])

 shear_mix.initNeighbourSearch("RotSphere",2.5,0.08)

 shear_mix.setNumTimeSteps(nt)
 shear_mix.setTimeStepSize(dt)
 shear_mix.readGeometry(geometry_name)

 # particle interactions
 #frictional interactions if they aren't bonded - elastic if they are
 fric_coeff=0.5
 stiffness=1.0
 fric_coeffw=0.1

 fric_prm=RotFrictionPrms("friction",stiffness,fric_coeff,fric_coeff,stiffness)
 fric_prmwall1=RotFrictionPrms("friction_wall1",stiffness,fric_coeffw,fric_coeffw,stiffness)
 fric_prmwall2=RotFrictionPrms("friction_wall2",stiffness,fric_coeffw,fric_coeffw,stiffness)
 fric_prm_marker=RotFrictionPrms("friction_marker",stiffness,fric_coeffw,fric_coeffw,stiffness)
 fric_prm_m_g=RotFrictionPrms("frict_mark_gouge",stiffness,fric_coeffw,fric_coeffw,stiffness)

# 30 is gouge
# 14 is a vertical line in the gouge to tract deformation
# 0 is the bonded particles for the wall

 shear_mix.createInteractionGroupTagged(fric_prm,30,-1,30,-1) #s-s
 shear_mix.createInteractionGroupTagged(fric_prm_marker,14,-1,14,-1) #s-s
 shear_mix.createInteractionGroupTagged(fric_prm_m_g,14,-1,30,-1)
 shear_mix.createInteractionGroupTagged(fric_prmwall1,14,-1,0,-1) #weak-wall
 shear_mix.createInteractionGroupTagged(fric_prmwall2,30,-1,0,-1) #wall gouge=weak

 # Make bonded interactions
 k=1.0
 kbreak=1000.0
 bond_param=RotBondPrms("bonded",
 k,
 k*0.2,
 k*0.08,
 k*0.035,
 kbreak*0.00025,
 kbreak*0.00125,
 kbreak*0.000125,
 kbreak*0.000125,
 0)
 shear_mix.createInteractionGroup(bond_param)

 shear_mix.createExclusion("bonded","friction")

 #create upper and lower walls to move and compress the model
 shear_mix.createWall(
 name="lowerWall",
 posn=Vec3(0.0,0.0,0.0),
 normal=Vec3(0.0,1.0,0.0)
 )
 shear_mix.createWall(
 name="upperWall",
 posn=Vec3(0.0,ymax,0.0),
 normal=Vec3(0.0,-1.0,0.0)
 )

 #side walls to keep the particles from falling out
 shear_mix.createWall(
 name="frontWall",
 posn=Vec3(0.0,0.0,0.0),
 normal=Vec3(0.0,0.0,1.0)
 )
 shear_mix.createWall(
 name="backWall",
 posn=Vec3(xmax,ymax,zmax),
 normal=Vec3(0.0,0.0,-1.0)
 )

 #Make walls interact with grains, first two elastic
 wp1=NRotElasticWallPrms("frontWallInteraction","frontWall",1.0)
 wp2=NRotElasticWallPrms("backWallInteraction","backWall",1.0)

 #Make walls bonded to grains
 wp3=NRotBondedWallPrms("upperWallInteraction","upperWall",1.0,4)
 wp4=NRotBondedWallPrms("lowerWallInteraction","lowerWall",1.0,8)

 #Create particle wall interactions
 shear_mix.createInteractionGroup(wp1)
 shear_mix.createInteractionGroup(wp2)
 shear_mix.createInteractionGroup(wp3)
 shear_mix.createInteractionGroup(wp4)

 shear_mix.createInteractionGroup(LocalDampingPrms(name="damping1",viscosity=0.6))
 shear_mix.createInteractionGroup(RotLocalDampingPrms(name="damping2",viscosity=0.6))

 kin_file_name="e_kin.dat"
 ek_prm=ParticleScalarFieldSaverPrms(
 fieldName="e_kin",
 fileName=kin_file_name,
 fileFormat="SUM",
 beginTimeStep=0,
 endTimeStep=nt,
 timeStepIncr=1
 )
 shear_mix.createFieldSaver(ek_prm)

 # in matlab add, line by line, epot, epot_w1, epot_w2.dat. Energy at each timestep
 ep_prm=InteractionScalarFieldSaverPrms(
    interactionName="friction",
    fieldName="potential_energy",
    fileName="epot.dat",
    fileFormat="SUM",
    beginTimeStep=0,
    endTimeStep=nt,
    timeStepIncr=1
 )
 shear_mix.createFieldSaver(ep_prm)

 # Number of frictional interactions for each kind of frictional particle
 # either add them line by line - or leave apart for analysis of each one
 nf_prm=InteractionScalarFieldSaverPrms(
 interactionName="friction",
 fieldName="count",
 fileName="nfriction.dat",
 fileFormat="SUM",
 beginTimeStep=0,
 endTimeStep=nt,
 timeStepIncr=1
 )
 shear_mix.createFieldSaver(nf_prm)

 nf_prm=InteractionScalarFieldSaverPrms("bonded","count","nbonds.dat","SUM",0,nt,1)
 shear_mix.createFieldSaver(nf_prm)

 wf_prm=WallVectorFieldSaverPrms(["lowerWall","upperWall"],"Force","wforce.dat","RAW_SERIES",0,nt,1)
 shear_mix.createFieldSaver(wf_prm)

 #wall position saver
 wp_prm=WallVectorFieldSaverPrms(["lowerWall","upperWall"],"Position","wposition.dat","RAW_SERIES",0,nt,1)
 shear_mix.createFieldSaver(wp_prm)

  # Create a checker to save snap shots
 cp=CheckPointPrms(
 fileNamePrefix="snaps/cpt", #send these to directory snaps
 beginTimeStep=0,
 endTimeStep=nt,
 timeStepIncr=dt_snap
 )
 #instantiante checkpointer
 shear_mix.createCheckPointer(cp)

 #output force chains for strong particles
 partforce_prm=InteractionVectorFieldSaverPrms(
    interactionName="friction",
    fieldName="force",
    fileName="forces/particleforce",
    fileFormat="RAW2",
    beginTimeStep=0,
    endTimeStep=nt,
    timeStepIncr=dt_force
    )
 shear_mix.createFieldSaver(partforce_prm)

 partforce_prm=InteractionVectorFieldSaverPrms(
    interactionName="frict_mark_gouge",
    fieldName="force",
    fileName="forces/particleforce_mg",
    fileFormat="RAW2",
    beginTimeStep=0,
    endTimeStep=nt,
    timeStepIncr=dt_force
    )
 shear_mix.createFieldSaver(partforce_prm)

 partforce_prm=InteractionVectorFieldSaverPrms(
    interactionName="friction_marker",
    fieldName="force",
    fileName="forces/particleforce_mark",
    fileFormat="RAW2",
    beginTimeStep=0,
    endTimeStep=nt,
    timeStepIncr=dt_force
    )
 shear_mix.createFieldSaver(partforce_prm)

 # Save the simulation to restart
 shear_mix.createRestartCheckPointer(
  CheckPointPrms(
  fileNamePrefix="saver/SAVE",
  beginTimeStep=dt_save,
  endTimeStep=nt,
  timeStepIncr=dt_save
  )
 )

 vrs=shear_mix.getLsmVersion()

 print "Current ESYS Version: ", vrs

 load=Loading(shear_mix,dy,t0,fn)
 shear_mix.addPreTimeStepRunnable(load)

 # execute the simulation
 shear_mix.run()

if (__name__ == "__main__"):
 runSimulation()

Hi Andy,

If you figure out the source of the irreproducibility, please post it here. It would be nice to know what causes nasty script-level problems like that.

Cheers,

Dion.

Hello all-

It seems when models diverge it is the result of a single timestep drop and recovery of the normal stress. Sometimes I have a a drop on the normal force on one of the walls. It lasts 1 timestep and then recovers back to the imposed value. The reproducibility problem seems to occur as a result of this. I had assumed this was a writing problem etc, but maybe it actually happens and as a result some grains rearrange? What could cause it? Has anyone else observed this? Sometimes it never happens, sometimes I get 1 or up to ~5. It seems random in displacement and time, and could occur on either of the servo walls.

AndyR