model new

model largestrain off
model config dynamic fluid plugin
model dynamic active off
model fluid active off

;-- Create Grid --
zone create brick size 1 1 32

;-- Define groups --
zone group 'M1' range position-z 29.0 32.0
zone group 'M2' range position-z 27.0 29.0
zone group 'M3' range position-z 19.0 27.0
zone group 'M4' range position-z 13.0 19.0
zone group 'AC'      range position-z  5.0 13.0
zone group 'AGS'     range position-z  0.0  5.0

;-- Define Ranges --
model range create 'x_neg' position-x 0.0
model range create 'x_pos' position-x 1.0
model range create 'y_neg' position-y 0.0
model range create 'y_pos' position-y 1.0
model range create 'z_bas' position-z 0.0
;
zone gridpoint group "XY00" range position-x 0 position-y 0
zone gridpoint group "XY01" range position-x 0 position-y 1
zone gridpoint group "XY11" range position-x 1 position-y 1
zone gridpoint group "XY10" range position-x 1 position-y 0

[global masterGps =  list(gp.list)(gp.isgroup(::gp.list, "XY11"))]
[global slave1Gps =  list(gp.list)(gp.isgroup(::gp.list, "XY00"))]
[global slave2Gps =  list(gp.list)(gp.isgroup(::gp.list, "XY01"))]
[global slave3Gps =  list(gp.list)(gp.isgroup(::gp.list, "XY10"))]

;Boundary conditions
zone face apply velocity-z 0.0   range named-range 'z_bas'
zone face apply velocity-x 0.0   range named-range 'x_neg' named-range 'z_bas' not
zone face apply velocity-x 0.0   range named-range 'x_pos' named-range 'z_bas' not
zone face apply velocity-y 0.0   range named-range 'y_neg' named-range 'z_bas' not
zone face apply velocity-y 0.0   range named-range 'y_pos' named-range 'z_bas' not

; Free-field Condition
fish define attach_gp
    loop local ii (2,list.size(MasterGps))
        local gp1 = MasterGps(ii)
        local gp2 = Slave1Gps(ii)
        local gp3 = Slave2Gps(ii)
        local gp4 = Slave3Gps(ii)
        local a1 = attach.create(gp2, gp1)
        local a2 = attach.create(gp3, gp1)
        local a3 = attach.create(gp4, gp1)
        attach.snap(a1) = 0
        attach.snap(a2) = 0
        attach.snap(a3) = 0
    endloop
end
[attach_gp]

;-- Set gravity --
[global grav = 9.81]
[global kw = 1.0]
[global uww = grav*kw]
model gravity (0, 0, [-grav])

;-- Fluid model and properties --
zone fluid cmodel assign isotropic
zone initialize fluid-density [kw]
zone gridpoint initialize fluid-modulus  0
zone gridpoint initialize fluid-tension -100
; permeability = normal perm * g*rho_w
zone fluid property permeability = [3.0e-3/uww] porosity 0.40 range group 'M1'
zone fluid property permeability = [3.0e-3/uww] porosity 0.40 range group 'M2' 
zone fluid property permeability = [3.0e-3/uww] porosity 0.40 range group 'M3' 
zone fluid property permeability = [3.0e-3/uww] porosity 0.40 range group 'M4' 
zone fluid property permeability = [3.0e-6/uww] porosity 0.40 range group 'AC'
zone fluid property permeability = [3.0e-4/uww] porosity 0.40 range group 'AGS'

;-- Pore Pressures --
zone water density 1.0
zone water plane normal 0.0,0.0,-1.0 origin 0.0,0.0,29.0
zone gridpoint fix pore-pressure 0.0 range position-z 29.0 32.0
zone gridpoint ini saturation 0.0 range position-z 29.0 32.0
zone gridpoint ini saturation 1.0 range position-z  0.0 29.0

;------------------------------------------------------------------------------------;
;  Solve for gravity
;------------------------------------------------------------------------------------;
zone cmodel assign elastic
;
zone property density [15.3/grav] shear 4.51e4 bulk 9.77e4 range group 'M1'
zone property density [15.3/grav] shear 5.66e4 bulk 1.23e5 range group 'M2' 
zone property density [15.3/grav] shear 8.63e4 bulk 1.87e5 range group 'M3' 
zone property density [15.3/grav] shear 8.63e4 bulk 1.87e5 range group 'M4' 
zone property density [16.1/grav] shear 6.61e4 bulk 1.43e5 range group 'AC' 
zone property density [15.6/grav] shear 1.19e5 bulk 2.59e5 range group 'AGS' 
;
zone initialize-stresses ratio 0.6
;
;Solve for gravity
model dynamic active off
model fluid active off
model mechanical active on
;
model solve mechanical ratio 1e-5
model save 'gravity1.sav'
;
; Assign dafalias model
zone cmodel assign dafalias range union group 'M2' group 'M3' group 'M4'
zone property density 2 evd0 0.8 Patm 100.0 range cmodel 'dafalias'
zone property G0 125.0 poisson 0.05 Mc 1.25 c 0.712 lambda 0.019 ec0 0.934 xi 0.7 range cmodel 'dafalias'
zone property mm 0.02 h0 7.05 ch 0.968 nb 1.25 A0 0.704 nd 2.1 zmax 2 cz 600.0 range cmodel 'dafalias'
fish operator ini_stress(modelname, zp)
    if zone.model(zp) = modelname
        local pp = zone.pp(zp)
        zone.prop(zp,'sxx-ini') = zone.stress.xx(zp) + pp
        zone.prop(zp,'syy-ini') = zone.stress.yy(zp) + pp
        zone.prop(zp,'szz-ini') = zone.stress.zz(zp) + pp
        zone.prop(zp,'sxy-ini') = zone.stress.xy(zp)
        zone.prop(zp,'sxz-ini') = zone.stress.xz(zp)
        zone.prop(zp,'syz-ini') = zone.stress.yz(zp)
    endif 
end
[ini_stress('dafalias', ::zone.list)] 
;
model solve mechanical ratio 1e-5
model save 'gravity2.sav'
;
zone gridpoint initialize displacement (0,0,0)
zone gridpoint initialize displacement-small (0,0,0)
zone gridpoint initialize velocity (0,0,0)
zone initialize state 0
;
;-- Configure for dynamic analysis --
model mechanical active on
model dynamic active on
model fluid active on
;
model dynamic time-total=0.0
model fluid   time-total=0.0
zone gridpoint initialize fluid-modulus 2.0e6 range position-x 0.0 29.0

;-- Dynamic Multistepping --
zone dynamic multi-step on

; Damping
zone dynamic damping hysteretic hardin 0.30 range group 'AC' 
zone dynamic damping hysteretic hardin 0.07 range group 'AGS'
zone dynamic damping hysteretic hardin 0.03 range group 'Masado1'

zone dynamic damping rayleigh 0.0149 3.47 

; rigid base
table '0' import 'Kobe_PI_32m_00_acc.txt' 
zone face apply acc-x 1.0 table '0' range named-range 'z_bas'
zone face apply-remove velocity-x range named-range 'x_neg' named-range 'z_bas' not
zone face apply-remove velocity-x range named-range 'x_pos' named-range 'z_bas' not

;Save pre-shaking model conditions
model save 'preEQ2.sav'
model dynamic timestep fix 2.0e-5

history delete
history interval 50
model history name='time' dynamic time-total
zone history name='pp_23' pore-pressure source gridpoint position 0 0 23.0
zone history name='acc-x_00m' acc-x         position 0 0 32.0
zone history name='acc-x_16m' acc-x         position 0 0 16.0
zone history name='acc-x_32m' acc-x         position 0 0  0.0
zone history name='vel-x_00m' vel-x         position 0 0 32.0
zone history name='vel-x_16m' vel-x         position 0 0 16.0
zone history name='vel-x_32m' vel-x         position 0 0  0.0
zone history name='dis-x_00m' dis-x         position 0 0 32.0
zone history name='dis-x_16m' dis-x         position 0 0 16.0
zone history name='dis-x_32m' dis-x         position 0 0  0.0
;
model solve dynamic time-total 21
;
model save '21s.sav'
;