Python script execution in PFC

Hello,

I am trying to simulate the direct shear test using PFC (Python scripting). When I plot the graph of vertical strain and shear strain, it is not starting from zero. Please go through the code and do the needful in this regard. Please find the python script for DST below.

new
set random 10001
domain extent -0.2 0.2
[ballFriction=0.0]
[wallFriction=0.0]
cmat default model linear method deform emod 1.0e8 kratio 2.5

def setup
 height=0.053
 height_=-0.053
 width=0.053
 width_=-0.053
 new_height_1=height+0.02
 new_height_2=-(height+0.02)
 new_width_1=width+0.02
 new_width_2=-(width+0.02)
 new_width_3=-(width+0.05)
 new_width_4=(width+0.05)
 fc=0.5
 a_damp=0.7
 p_density=1200
end
@setup

def make_walls
 command
  wall create vertices @new_width_2, @height @new_width_1, @height id 1, name top_wall
  wall create vertices @new_width_2, @height_ @new_width_1, @height_ id 2, name bot_wall
  wall create vertices @width_, 0, @width_, @new_height_1 id 3
  wall create vertices @width, 0, @width, @new_height_1 id 4
  wall create vertices @width_, @new_height_2 @width_, 0 id 5
  wall create vertices @width, @new_height_2 @width, 0 id 6
  wall create vertices @new_width_3, 0 @width_, 0 id 7
  wall create vertices @width, 0 @new_width_4, 0 id 8
 end_command
end
@make_walls

wall property kn=1e8

def wll_wp
      global wpy1 = wall.find(2)
      global wpy2 = wall.find(1)
      global wpx1 = wall.find(3)
      global wpx2 = wall.find(4)
      global wpx11 = wall.find(5)
      global wpx22 = wall.find(6)
end
@wll_wp

ball distribute porosity 0.175 radius 0.001 0.004 box -0.053 0.053 -0.053 0.053
ball attribute density 1700 damp 0.7
ball property fric @ballFriction
wall property fric @wallFriction

cycle 10 calm 10
set timestep scale
solve arat 1e-5
set timestep auto
cycle 10
solve arat 1e-5

define compute_vessel_dimensions
  global et3_xlen = wall.pos.x(wall.find(4)) - wall.pos.x(wall.find(3))
  global et3_ylen = wall.pos.y(wall.find(1)) - wall.pos.y(wall.find(2))
  x_wpx1= wall.pos.x(wpx1)
  x_wpx2= wall.pos.x(wpx2)
  y_wpy1= wall.pos.y(wpy1)
  y_wpy2= wall.pos.y(wpy2)
end

def _et3_wstresses
         et3_wss
          _xfob = 0.5 * (wall.force.contact(wpx1,1)+wall.force.contact(wpx11,1) - wall.force.contact(wpx2,1)-wall.force.contact(wpx22,1))
          _yfob = 0.5 * (wall.force.contact(wpy1,2) - wall.force.contact(wpy2,2))
        if _et3_wax # 0.0
         et3_wsxx = _xfob/_et3_wax
        else
         et3_wsxx = 0.0
        end_if
        if _et3_way # 0.0
         et3_wsyy = _yfob/_et3_way
        else
         et3_wsyy = 0.0
        end_if
end

define et3_wss
    compute_vessel_dimensions
    _et3_way = et3_xlen*1.0
    _et3_wax = et3_ylen*1.0
end

define et3_servo_gain_poly
  if et3_servo_xon = 1 then
     sum_knx = 0.0
    clist = wall.contactmap(wpx1)
    loop foreach local cp clist
      sum_knx = sum_knx + contact.prop(cp,'kn')
    endloop
    clist = wall.contactmap(wpx11)
    loop foreach cp clist
      sum_knx = sum_knx + contact.prop(cp,'kn')
    endloop
    clist = wall.contactmap(wpx2)
    loop foreach cp clist
      sum_knx = sum_knx + contact.prop(cp,'kn')
    endloop
    clist = wall.contactmap(wpx22)
    loop foreach cp clist
      sum_knx = sum_knx + contact.prop(cp,'kn')
    endloop
    sum_knx = 0.5 * sum_knx
    if sum_knx # 0.0 then
      et3_servo_gx = et3_servo_alpha * _et3_wax / (sum_knx * global.timestep)
    else
      et3_servo_gx = 0.0
    end_if
  end_if
  ;
  if et3_servo_yon = 1 then
    sum_kny = 0.0
    clist = wall.contactmap(wpy1)
    loop foreach cp clist
      sum_kny = sum_kny + contact.prop(cp,'kn')
    end_loop
    clist = wall.contactmap(wpy2)
    loop foreach cp clist
      sum_kny = sum_kny + contact.prop(cp,'kn')
    end_loop
    sum_kny = 0.5 * sum_kny
    if sum_kny # 0.0 then
      et3_servo_gy = et3_servo_alpha * _et3_way / (sum_kny * global.timestep)
    else
      et3_servo_gy = 0.0
    end_if
  end_if
end

define et3_servo_poly
  _et3_wstresses
  if et3_servo_xon = 1
    local _sgn = math.sgn( et3_wsxx - et3_wsxx_req )
    if et3_wsxx # 0.0 then
      wall.vel(wpx1,1) = et3_servo_gx * (et3_wsxx - et3_wsxx_req)
      wall.vel(wpx11,1) = et3_servo_gx * (et3_wsxx - et3_wsxx_req)
    else
      wall.vel(wpx1,1) = _sgn * _vmax
      wall.vel(wpx11,1) = _sgn * _vmax
    end_if
    if math.abs(wall.vel(wpx1,1)) > _vmax then
      wall.vel(wpx1,1) = _sgn * _vmax
      wall.vel(wpx11,1) = _sgn * _vmax
    end_if
    wall.vel(wpx2,1) = -wall.vel(wpx1,1)
    wall.vel(wpx22,1) = -wall.vel(wpx1,1)
  end_if
  ;
  if et3_servo_yon = 1
    _sgn = math.sgn( et3_wsyy - et3_wsyy_req )
    if et3_wsyy # 0.0 then
      wall.vel(wpy1,2) = et3_servo_gy * (et3_wsyy - et3_wsyy_req)
    else
      wall.vel(wpy1,2) = _sgn * _vmax
    end_if
    if math.abs(wall.vel(wpy1,2)) > _vmax then
      wall.vel(wpy1,2) = _sgn * _vmax
    end_if
    wall.vel(wpy2,2) = -wall.vel(wpy1,2)
  end_if
  ;
end

define et3_servo_gain

  _et3_gain_cnt = _et3_gain_cnt + 1
  if _et3_gain_cnt >= et3_servo_gain_cyc then
    _et3_gain_cnt = 0
    et3_servo_gain_poly
  end_if
end

def m_tol
 y_tol=(et3_wsyy - et3_wsyy_req)/et3_wsyy_req
 x_tol=(et3_wsxx - et3_wsxx_req)/et3_wsxx_req
 mech_solve=mech.solve("aratio")
end

history id 104 @y_tol
history id 105 @x_tol
history id 107 @mech_solve

set fish callback 1.0 @m_tol

[et3_wsxx_req = -75.0e3]
[et3_wsyy_req = -75.0e3]
[et3_servo_xon=1.0]
[et3_servo_yon=1.0]

set fish callback 1.0 @et3_servo_poly

history id 101 @et3_wsxx
history id 102 @et3_wsyy
history id 103 @mech_solve

set orientation on
calm

[tol = 5e-3]
[_et3_gain_cnt = 0]
[et3_servo_gain_cyc = 50]
[et3_servo_alpha = 0.5]
[et3_servo_vmax = 10.0]
[_vmax=et3_servo_vmax]

define stop_me
  et3_servo_gain
  if math.abs((et3_wsyy - et3_wsyy_req)/et3_wsyy_req) > tol
    exit
  endif
  if math.abs((et3_wsxx - et3_wsxx_req)/et3_wsxx_req) > tol
    exit
  endif
  if mech.solve("aratio") > 1e-5
    exit
  endif
  command
   sav sample_ini.p2sav
  end_command
  stop_me = 1
end

set grav 0 0

solve fishhalt @stop_me

res sample_ini.p2sav

def ini_shear
 et3_xlen_0=et3_xlen
 et3_ylen_0=et3_xlen
end
@ini_shear

def get_ss
wall_x_dis=(wall.disp.x(wpy2)+wall.disp.x(wpx1)+wall.disp.x(wpx2))/3.0
wall_x_str=wall_x_dis/et3_ylen_0
wall_y_dis=(wall.disp.y(wpy2)+wall.disp.y(wpx1)+wall.disp.y(wpx2))/3.0
wall_y_str=wall_y_dis/et3_ylen_0
Fx1=wall.force.contact(wpx1,1)
Fx2=wall.force.contact(wpx2,1)
Fx=Fx1+Fx2
stress_x = -Fx/et3_xlen_0
vol_str = (et3_xlen*et3_ylen-et3_xlen_0*et3_ylen_0)/(et3_xlen_0*et3_ylen_0)
end

define stop_shear
  gain_cnt = gain_cnt + 1
  if gain_cnt >= et3_servo_gain_cyc then
    et3_servo_gain_poly
    gain_cnt = 0
  endif
 if math.abs(wall_x_str) >= target then
    stop_shear = 1
  endif
end

set fish callback 1.0 @get_ss

history id 1 @wall_x_dis
history id 2 @wall_x_str
history id 3 @stress_x
history id 4 @vol_str
history id 5 @wall_y_str

[et3_servo_xon=0]

ball prop fric 0.5
[rate = 1.0e-3]
[target = 2.0e-1]
[stop_shear = 0]
[gain_cnt=0]

wall attribute xvelocity [rate] range id 1
wall attribute xvelocity [rate] range id 3
wall attribute xvelocity [rate] range id 4
wall attribute xvelocity [rate] range id 7

set echo off
solve fishhalt @stop_shear
set echo on