; compl12 is similar to compl11 but improved resolution to 300*200. 
; The permeability contrast between seal and permeable layer has been decreased
; by increasing seal permeability to 1e-16 (from 1e-18).
; A further difference is that when both material yields the increase in permeability
; is less intense than in compl11 (1e-14, 5e-13, 1e-13 for the different states)

; compl8 equal compl5 intermediate velocity but shear and bulk have been further increased in the imperm. layer. and
; decreased p at top to 1 kbar
; compl9  kept the same moduli and Mohr-Coulomb parameters, reduced the veloc by half and added a permability function from
; Nerys's pump, 

tit
FLAC Modelling Input File
;
config extra 4 gw
grid  300   200 
gen 0.00e0,0.00e0 0.00e0,1.00e4 2.8e4,1.00e4 2.8e4,0.0e0 i=1,301 j=1,201
mark i=1
mark j=1
mark i = 301 
mark j = 201


***************************************
model mohr
water bulk 2e9 dens 1000
prop poro=0.3 ; to be elaborated below
;; 1 mech step = 1 fluid step; modify if necessary!!!
set nmech 1 ngw 10 flow on
***************************************

tab 1 0, 0, 0, 10000, 12500, 10000, 15500, 5000, 28000, 5000, 28000, 0

gen tab 1


;
; Properties for each sub-regions: 
pro dens= 2720  bulk=4.0e9 shear=2.45e9 reg=20, 20
pro coh=3.00e7 tens=1.5e7 reg=20, 20
pro fric=40 dil= 10  reg=20, 20                     ; changed fric and dil compared to all previour "compl" runs
pro perm=1.00e-14 porosity=0.2 reg=20, 20

; Properties for each sub-regions: 
pro dens= 2720  bulk=8.e10 shear=2.45e10 reg=250, 150 ; bulk and shear are further increased compared to compl5 
pro coh=3.00e7 tens=1.5e7 reg=250, 150
pro fric=40 dil= 10  reg=250, 150
pro perm=1.00e-16 porosity=0.2 reg=250, 150

plo hold perm blo

;
apply syy = -1.e8 j= 201 
*****
ini sxx=-3.66800e8 var=0, 2.668e8 i=1, 301 j=1,201
ini syy=-3.66800e8 var=0, 2.668e8 i=1, 301 j=1,201
ini szz=-3.66800e8 var=0, 2.668e8 i=1, 301 j=1,201
ini sxy=0e6


fix x i=1
fix x i= 301 
fix y j=1

;
def hydro_pp
loop i (1,igp)
loop j (1,jgp)
if j = jgp Then
sss2 = 3.67e7
;sss2 = 7.35e7
else
aaa = abs(y(i, j) - y(i, jgp))
sss2 = 1000 * 9.8 * aaa+3.67e7
;sss2 = 1000 * 9.8 * aaa+7.35 e7 ;for 2kb
end_if
gpp(i, j) = sss2
end_loop
end_loop
end
hydro_pp
*
;from Nerys's pump
def store_perm
    loop i (1,izones)
      loop j (1,jzones)
        ex_1(i,j) = k11(i,j)
      end_loop
    end_loop
end
store_perm
*

def perm_yield
;take this out for one test
while_stepping
    
    loop i (1,300)
        loop j (1,200)
               if state(i,j)=1 then ; at yield in shear or volume
                
                  k11(i,j) = 5.e-13
                  k12(i,j) = 0.0
                  k22(i,j) = 5.e-13
                
              else
              if state(i,j)=2 then ; yielded in the past, not now
               
                  k11(i,j) = 1.e-13
                  k12(i,j) = 0.0
                  k22(i,j) = 1.e-13
               
              else
              if state(i,j)=3 then ; yield in tension
               
                  k11(i,j) = 5.e-13
                  k12(i,j) = 0.0
                  k22(i,j) = 5.e-13
               
              end_if
              end_if
              end_if
              end_loop
    end_loop
end
perm_yield
;
ini pp=3.67e7 j=201
fix pp j=201

ini sat 1
fix sat
;
set gravity=9.8
;
hist 1 unbal nstep 1

set small
set mech on fl off
step 200
ini   xd 0    yd 0
ini xvel 0  yvel 0
;plot hold hist 1 
set mech on fl on
step 400  ; this can be increased if you think the unbalanced forces haven't bottomed out
ini   xd 0    yd 0 
ini xvel 0  yvel 0

set large
plot hold hist 1


;initial velocity command should be here ie ini xvel 0.004 etc;


;stress loading time
set nmech 1 ngw 1 flow on
ini xvel  0.001 var=-0.002 ,0 i=1, 301  j=1, 201 
step 5000
save compl12_0001.sav
step 5000
save compl12_0002.sav



set nmech 1 ngw 10 flow on


;ini xvel  0.001 var=-0.002 ,0 i=1, 301  j=1, 201 ; halved the velocity compared to compl8

;set mech off
apply dis= 1e-10 i= 1 301 j= 1


def files
    if nn<10 then
     savefile = 'compl12_00' + string(nn) + '.sav'
    else
    if nn>9 then
     if nn<100 then
     savefile = 'compl12_0' + string(nn) + '.sav'
    else
    if nn>99 then
     savefile = 'compl12_' + string(nn) + '.sav'
    end_if
    end_if
    end_if
    end_if
    COMMAND
        save @savefile
    END_COMMAND
end

define step_loop
loop nn (1,60)
   ;loop n (1,1000)
	command
 		step 10000
	end_command
      ;vsi_perm
   ;end_loop
files
end_loop
end
step_loop




return