;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;The code is to remove the right wall instantaneously to trigger the collapse of the initial granular column with ; 
;fractal dimension 3.5 (IGC3p5). It is a main code and other functions can be added based on it.                  ;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

New

restor IGC3p5.sav

set disk on 

title IGCFlow3p5

dele wal 2

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;The right wall is then removed. The simulation parameters are specified;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

def fric_parameter
 wbfric = 0.60
end
fric_parameter

wal id 2 nodes (L_1,0)   (1000,0)
wall prop  id 1 kn bkn ks bks  fric wbfric
wall prop  id 2 kn bkn ks bks  fric wbfric
prop dens bdens kn bkn ks bks  fric wbfric

damp local 0.0
damp viscous normal 0.12
damp viscous shear  0.12

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;The small, medium and large sized particles are described by three different colors, respectively;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

def colorball

  bp = ball_head
  loop while bp # null
   
    brr = b_rad(bp)

    if brr < 0.40 then
     b_color(bp) = 0
    endif

    if brr = 0.40 then
     b_color(bp) = 0
    endif

    if brr = 0.60 then
     b_color(bp) = 1
    endif

    if brr = 0.80 then
     b_color(bp) = 1
    endif

    if brr = 1.00 then
     b_color(bp) = 2
    endif

    if brr = 1.20 then
     b_color(bp) = 2
    endif

  bp = b_next(bp)
  end_loop

end
colorball

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;The physical motion time of granular flows is set and recorded;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

def setTime
 oldtime = time
end
setTime

def getTime
 getTime = time - oldTime
end
getTime

hist gettime  ;1 hist number

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;The simulation files at the normalized time of 1.0, 2.0, 3.0, 4.0, 5.0; 
;6.0 7.0 and 8.0 are specified physical time are saved, respectively.  ;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

def pan_duan

 WHILESTEPPING

 pdd_1 = abs(gettime - 2.6133)   ;The normalized time is 1.0
 pdd_2 = abs(gettime - 5.2267)   ;The normalized time is 2.0
 pdd_3 = abs(gettime - 7.8401)   ;The normalized time is 3.0
 pdd_4 = abs(gettime - 10.4535)  ;The normalized time is 4.0
 pdd_5 = abs(gettime - 13.0665)  ;The normalized time is 5.0
 pdd_6 = abs(gettime - 15.6798)  ;The normalized time is 6.0
 pdd_7 = abs(gettime - 18.2931)  ;The normalized time is 7.0
 pdd_8 = abs(gettime - 20.9064)  ;The normalized time is 8.0

  if  pdd_1 < 0.0001
    command
      sav 3p5_normalized_1.sav
    endcommand
  end_if 

  if  pdd_2 < 0.0001
    command
      sav 3p5_normalized_2.sav
    endcommand
  end_if 

  if  pdd_3 < 0.0001
    command
      sav 3p5_normalized_3.sav
    endcommand
  end_if 

  if  pdd_4 < 0.0001
    command
      sav 3p5_normalized_4.sav
    endcommand
  end_if 

  if  pdd_5 < 0.0001
    command
      sav 3p5_normalized_5.sav
    endcommand
  end_if 

  if  pdd_6 < 0.0001
    command
      sav 3p5_normalized_6.sav
    endcommand
  end_if 

  if  pdd_7 < 0.0001
    command
      sav 3p5_normalized_7.sav
    endcommand
  end_if 

  if  pdd_8 < 0.0001
    command
      sav 3p5_normalized_8.sav
    endcommand
  end_if 
end
pan_duan

cy 99999999999999999