; filename: make_brick.dat
; 
; Generate an isotropic, homogeneous assembly of balls with periodic
; boundary conditions, and settle under low isotropic confinement 
; using a servo-control algorithm to distort the periodic domain 
;-------------------------------------------------------------------------initial parameter
model new
; Generate packing
model random 10001
[resolution = 23]   ; approx. number of balls accross specimen width 
[davg  = 2.5e-2]    ; average ball diameter (m)
[dratio = 1.66]     ; ball max/min diameter ratio [-]
[tag = 'A'+ string(resolution)]

;-------------------------------------------------------------------------create domain
; Setup domain extent (aspect ratio 2:1) and boundary conditions
model domain extent [-0.5*resolution*davg] [0.5*resolution*davg] ...
                    [-0.5*resolution*davg] [0.5*resolution*davg] ...
                    [-1.0*resolution*davg] [1.0*resolution*davg] ...
                     condition periodic
                     
; Distribute balls to a target low porosity and set density and 
; local damping attributes
[dmin = 2.0*davg / (1.0 + dratio)]
[dmax = dmin * dratio]
ball distribute porosity 0.5 radius [0.5*dmin] [0.5*dmax] ;loose sample
ball attribute density 2000.0 damp 0.7
;-------------------------------------------------------------------------contact model
contact cmat default model linear ...
                     method deformability emod 7e8 kratio 5.0 property fric 0.5 dp_nratio 0.22 dp_sratio 0.22
cmat apply
;Monitoring
model history name 'ratio-average' mechanical ratio-average
model history name 'timestep' timestep
model cycle 1000 calm 100
model solve ratio-average 5e-3
;--------------------------------------------------------------------------relax model
; Load utility functions
program call 'servo-triaxial.fis' suppress

; Parameter definitions (used by the servo-control function)
[servo_do_z = true]
[servo_do_x = true]
[servo_do_y = true]
[servo_ts_xx = -1e5  ]
[servo_ts_yy = -1e5  ]
[servo_ts_zz = -1e5  ]
[servo_gain  =  1.0e-5 ]
[servo_tol = 1e-2]
[s_xx = 0.0]
[s_yy = 0.0]
[s_zz = 0.0]

[Lx = domain.max.x() - domain.min.x()]
[Ly = domain.max.y() - domain.min.y()]
[Lz = domain.max.z() - domain.min.z()]
measure create id 1 radius [0.95*0.5*Ly]

;----------------------------------------------
; Monitoring
fish history name 's_xx' @s_xx
fish history name 's_yy' @s_yy
fish history name 's_zz' @s_zz
fish history name 'L_x'  @Lx
fish history name 'L_y'  @Ly
fish history name 'L_z'  @Lz

model domain tolerance 1e-3
model solve fish-call 1.0 @servo_domain ...
            fish-halt @servo_halt ...
            ratio-average 1e-4 ...
            cycles 100 and

model save "compaction_sample.sav"

program return
; EOF: make_brick.dat