#===============================================================================================================
#  Undrained Triaxial Compression Test
#
#  REFERENCE PAPER: CFD-DEM model to assess stress-induced anisotropy in undrained granular material
#  Latest Update: 08-DEC-2020
#  
#  PROGRESS
#	1. Sample preparation
#		1.1 Sample dimension:	0.02 x  0.02 x 0.04 
#		1.2 No. of particle:	target->18000	obtained-> 7073
#		1.3 Void ratio: 	target->0.65	obtained-> 0.584
#		1.4 Particle dist. maintained as per ref.:	yes
#		1.5 Simulation time for sample filling:		In code->	Actual->
#	2. Confinment
#		2.1 sigma_3: 160 MPa
#		2.2 Inertial Factor:
#
#===============================================================================================================

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#	Common Setting for restart file
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
atom_style	granular
atom_modify	map array
boundary	f f f
newton		off
#hard_particles 	yes
#soft_particles	yes
communicate	single vel yes
processors	* * *
units		si

region		reg block 0.0 2.8 0.0 2.8 0.0 4.8 units box	
region		reg_ins block 0.4 2.4 0.4 2.4 0.4 4.4 units box
	
#create_box	2 reg							# Comment it while reading restart file	
read_restart	restart.sample_settled

neighbor	0.01 bin
neigh_modify	delay 0 page 1000000 
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

#echo screen 
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#      Global Input Data
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Bulk_material: Sand
variable solid_density equal 1692	# kg/m3
variable sample_width equal 2.0		# m
variable sample_height equal 4.0	# m

# interaction property [NOTE: Stiffness value is same for both particle and wall]
# Bulk_material: Sand
variable kn_pp equal 6.24e7		# Normal stifness (N/m) particle-particle
variable kt_pp equal 2.496e6		# Tangential stifness (N/m) particle-particle
variable gn_pp equal 0.0		# Normal Damping coefficient particle-particle
variable gt_pp equal 0.0		# Tangential Damping coefficient particle-particle

# Interaction property Particle-wall
variable kn_pw equal 6.24e7		# Normal stifness (N/m) particle-wall
variable kt_pw equal 0.0		# Tangential stifness (N/m) particle-wall
variable gn_pw equal 0.0		# Normal Damping coefficient particle-wall
variable gt_pw equal 0.0		# Tangential Damping coefficient particle-wall

# interaction property Particle-particle 
variable us_pp equal 0.577		# Coefficient of static friction particle-particle
variable ur_pp equal 0.10		# Coefficient of static rolling particle-particle
variable uvd_pp equal 1.0		# Coefficient of static rolling Viscous damping particle-particle
variable beta equal 5.0			# Coefficient of static rolling stifness particle-particle

# interaction property Particle-wall
variable us_pw equal 0.0		# Coefficient of static friction particle-wall
variable ur_pw equal 0.0		# Coefficient of static rolling particle-wall
variable uvd_pw equal 0.0		# Coefficient of static rolling Viscous damping particle-particle

# triaxial test
variable time_step equal 1e-6
variable nstep equal 1000		# interval in which files are saved

variable vel  equal 500.0		# Maximum velocity
variable kp equal 0.0001	 	# proportional constant for PID controller
variable sigma_3 equal 100000		# Confining pressure [Pa]
variable sigma_1_min equal 110000	# Confining pressure [Pa] 
variable sigma_1_max equal 200000	# Confining pressure [Pa]
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

timestep      	${time_step}
#fix           	gravi all gravity 9.81 vector 0.0 0.0 -1.0

#shell rm -r post graph
shell mkdir post graph restart

# Particle distributions
# 15485863, 15485867, 32452843, 32452867, 49979687, 49979693, 67867967, 67867979, 86028121, 86028157, 98653427, 96835421
fix  pts1 all particletemplate/sphere 15485863 atom_type 1 density constant ${solid_density} radius constant 0.0407 #0.045
fix  pts2 all particletemplate/sphere 15485867 atom_type 1 density constant ${solid_density} radius constant 0.0416 #0.05
fix  pts3 all particletemplate/sphere 32452843 atom_type 1 density constant ${solid_density} radius constant 0.043 #0.055
fix  pts4 all particletemplate/sphere 32452867 atom_type 1 density constant ${solid_density} radius constant 0.0453 #0.06
fix  pts5 all particletemplate/sphere 49979687 atom_type 1 density constant ${solid_density} radius constant 0.0498 #0.075 
fix  pts6 all particletemplate/sphere 49979693 atom_type 1 density constant ${solid_density} radius constant 0.0543 #0.08
fix  pts7 all particletemplate/sphere 67867967 atom_type 1 density constant ${solid_density} radius constant 0.0679 #0.09

fix  pdd1 all particledistribution/discrete 86028157  7 pts1 0.28 pts2 0.25 pts3 0.18 pts4 0.12 pts5 0.08 pts6 0.05 pts7 0.03 

# Material properties required for new pair styles
fix 	p1 all property/global kn peratomtypepair 2 ${kn_pp} ${kn_pw} ${kn_pw} ${kn_pw}
fix 	p2 all property/global kt peratomtypepair 2 ${kt_pp} ${kt_pw} ${kt_pw} ${kt_pw} 
fix 	p3 all property/global gamman peratomtypepair 2 ${gn_pp} ${gn_pw} ${gn_pw} ${gn_pw}   
fix 	p4 all property/global gammat peratomtypepair 2 ${gt_pp} ${gt_pw} ${gt_pw} ${gt_pw}
fix	p5 all property/global coefficientFriction peratomtypepair 2 ${us_pp} ${us_pw} ${us_pw} ${us_pp}
fix 	p6 all property/global coefficientRollingFriction peratomtypepair 2 ${ur_pp} ${ur_pw} ${ur_pw} ${ur_pp}  
fix 	p7 all property/global coefficientRollingViscousDamping peratomtypepair 2 ${uvd_pp} ${uvd_pw} ${uvd_pw} ${uvd_pp}
fix     p8 all property/global rollingStiffness scalar ${beta}
fix 	p9 all property/global k_finnie peratomtypepair 2 1.0 1.0 1.0 1.0


# Importing CAD model-------------------------------------------------------------------------------------------
fix	boxwalls_x1 all wall/gran model hooke/stiffness tangential history rolling_friction epsd3 primitive type 2 xplane 2.4	
fix	boxwalls_x2 all wall/gran model hooke/stiffness tangential history rolling_friction epsd3 primitive type 2 xplane 0.4
fix	boxwalls_y1 all wall/gran model hooke/stiffness tangential history rolling_friction epsd3 primitive type 2 yplane 2.4
fix	boxwalls_y2 all wall/gran model hooke/stiffness tangential history rolling_friction epsd3 primitive type 2 yplane 0.4
fix	boxwalls_z1 all wall/gran model hooke/stiffness tangential history rolling_friction epsd3 primitive type 2 zplane 0.4
fix	boxwalls_z2 all wall/gran model hooke/stiffness tangential history rolling_friction epsd3 primitive type 2 zplane 4.4


# PARTICLE INSERTION--------------------------------------------------------------------------------------------
fix    ins all insert/pack seed 98653427  distributiontemplate pdd1 &
       maxattempt 100 insert_every once overlapcheck yes all_in yes vel constant 0.0 0.0 0.0 &
       region reg_ins particles_in_region 10648

# ----------------------------------CONTACT MODEL----------------------------------------------------------------
# For particle-particle interaction
pair_style  	gran model hooke/stiffness tangential history rolling_friction epsd3 
pair_coeff 	 * *
#----------------------------------------------------------------------------------------------------------------

fix	      integr all nve/sphere
compute       1 all erotate/sphere
compute       ke all ke/atom
compute       ke1 all ke
compute       overlap all contact/atom
compute	      fc all pair/gran/local pos force

thermo_style  custom time step atoms ke c_1 vol
thermo        ${nstep}
thermo_modify lost ignore norm no

fix  	energy all ave/time ${nstep} 1 ${nstep} c_ke1 c_1 file graph/energy.dat 

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#         Particle Insertion
#
#	(1) Insert particle by gravity filling
#	(2) Relax the sample till KE reduces to zero
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
run   	1
#dump  	dmp all custom/vtk ${nstep} post/particle_*.vtk id type type x y z ix iy iz vx vy vz fx fy fz omegax omegay omegaz radius c_ke c_overlap
#dump	dmp_fc all local/gran/vtk ${nstep} post/forcechain_*.vtk fc
unfix  	ins


#---------------------------------------------------------------------------
#	Radius Expansion method to grow the particle 
#---------------------------------------------------------------------------
compute         rad all property/atom radius
variable        dgrown atom 1.01*2.0*c_rad
fix             grow all adapt 10000 atom diameter v_dgrown
#run     	320000 						# Comment it while reading restart file
unfix		grow

#---------------------------------------------------------------------------
#	Relax Sample
#---------------------------------------------------------------------------

#run	200000

#write_restart restart.sample_settled				# Comment it while reading restart file	

#unfix	gravi	
unfix	boxwalls_x1 
unfix	boxwalls_x2 
unfix	boxwalls_y1 
unfix	boxwalls_y2 
unfix	boxwalls_z1
unfix	boxwalls_z2


region		del_region2 block 0.4 2.4 0.4 2.4 0.4 4.4  units box side out	
delete_atoms region del_region2

#===============================================================================================================

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#         Consolidation
#
#	(1) Remove the old wall and import new wall in the same position
#	(2) The new wall are made of servo type so that constant normal forces can ne maintained
#	(3) Change the wall/grain model as previous wall/grain command has old mesh detail  
#	(4) While using unfix and undump command maintain the bottom to top sequence
#	(5) Remove all graph data as it may get lost
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

variable f_top equal ${sigma_3}*${sample_width}*${sample_width}
variable f_side equal ${sigma_3}*${sample_width}*${sample_height}
variable kp1 equal 2*${kp}
variable nstep1 equal 1

fix	twall_servo all mesh/surface/stress/servo file meshes/topwall.stl type 2 com 1.4 1.4 4.4 ctrlPV force axis 0.0 0.0 1.0 target_val -${f_top} vel_max ${vel} kp ${kp1}
fix	bwall_servo all mesh/surface/stress/servo file meshes/bottomwall.stl type 2 com 1.4  1.4  0.4 ctrlPV force axis 0.0 0.0 1.0 target_val ${f_top} vel_max ${vel} kp ${kp1}
fix	wall1_servo all mesh/surface/stress/servo file meshes/sidewall1.stl type 2 com 2.4  1.4  2.4 ctrlPV force axis 1.0 0.0 0.0 target_val -${f_side} vel_max ${vel} kp ${kp}
fix	wall2_servo all mesh/surface/stress/servo file meshes/sidewall2.stl type 2 com 1.4  0.4  2.4 ctrlPV force axis 0.0 1.0 0.0 target_val ${f_side} vel_max ${vel} kp ${kp}
fix	wall3_servo all mesh/surface/stress/servo file meshes/sidewall3.stl type 2 com 0.4  1.4  2.4 ctrlPV force axis 1.0 0.0 0.0 target_val ${f_side} vel_max ${vel} kp ${kp}
fix	wall4_servo all mesh/surface/stress/servo file meshes/sidewall4.stl type 2 com 1.4  2.4  2.4 ctrlPV force axis 0.0 1.0 0.0 target_val -${f_side} vel_max ${vel} kp ${kp}

fix 	granwall2 all wall/gran model hooke/stiffness tangential history rolling_friction epsd3 mesh n_meshes 6 meshes twall_servo bwall_servo wall1_servo wall2_servo wall3_servo wall4_servo

fix  	force1 all ave/time ${nstep1} 1 ${nstep1} f_twall_servo[1] f_twall_servo[2] f_twall_servo[3] f_twall_servo[7] f_twall_servo[8] f_twall_servo[9] file graph/force_twall.dat 
fix  	force2 all ave/time ${nstep1} 1 ${nstep1} f_bwall_servo[1] f_bwall_servo[2] f_bwall_servo[3] f_bwall_servo[7] f_bwall_servo[8] f_bwall_servo[9] file graph/force_bwall.dat
fix  	force3 all ave/time ${nstep1} 1 ${nstep1} f_wall1_servo[1] f_wall1_servo[2] f_wall1_servo[3] f_wall1_servo[7] f_wall1_servo[8] f_wall1_servo[9] file graph/force_wall1.dat
fix  	force4 all ave/time ${nstep1} 1 ${nstep1} f_wall2_servo[1] f_wall2_servo[2] f_wall2_servo[3] f_wall2_servo[7] f_wall2_servo[8] f_wall2_servo[9] file graph/force_wall2.dat
fix  	force5 all ave/time ${nstep1} 1 ${nstep1} f_wall3_servo[1] f_wall3_servo[2] f_wall3_servo[3] f_wall3_servo[7] f_wall3_servo[8] f_wall3_servo[9] file graph/force_wall3.dat
fix  	force6 all ave/time ${nstep1} 1 ${nstep1} f_wall4_servo[1] f_wall4_servo[2] f_wall4_servo[3] f_wall4_servo[7] f_wall4_servo[8] f_wall4_servo[9] file graph/force_wall4.dat

#dump	dump1 all mesh/gran/VTK ${nstep} post/mesh_top_*.vtk stress twall_servo 
#dump	dump2 all mesh/gran/VTK ${nstep} post/mesh_bottom_*.vtk stress  bwall_servo 
#dump	dump3 all mesh/gran/VTK ${nstep} post/mesh_side1_*.vtk stress wall1_servo 
#dump	dump4 all mesh/gran/VTK ${nstep} post/mesh_side2_*.vtk stress wall2_servo 
#dump	dump5 all mesh/gran/VTK ${nstep} post/mesh_side3_*.vtk stress wall3_servo 
#dump	dump6 all mesh/gran/VTK ${nstep} post/mesh_side4_*.vtk stress wall4_servo

#restart  ${nstep} restart/restart_confinment_*.file

run	1000000				# Comment it while reading restart file


#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# 	Appling Deviator Stress [Resilant Modulus]
#	(1) modify the target value for the fix servo wall
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

variable load_step equal 0.05/${time_step}
variable unload_step equal 0.01/${time_step}
variable relax_step equal 0.9/${time_step}
variable kp_t equal 0.0020
variable kp_t1 equal 0.0001
variable kp_s equal 0.01

variable f_max equal ${sigma_1_max}*${sample_width}*${sample_width}
variable f_min1 equal 110000*${sample_width}*${sample_width}
variable f_min2 equal 120000*${sample_width}*${sample_width}
variable f_min3 equal 140000*${sample_width}*${sample_width}
variable f_min4 equal 160000*${sample_width}*${sample_width}
variable f_min5 equal 180000*${sample_width}*${sample_width}

fix_modify twall_servo servo/ctrlParam  ${kp_t} 0.0 0.0
fix_modify bwall_servo servo/ctrlParam  ${kp_t} 0.0 0.0

fix_modify wall1_servo servo/ctrlParam  ${kp_s} 0.0 0.0
fix_modify wall2_servo servo/ctrlParam  ${kp_s} 0.0 0.0
fix_modify wall3_servo servo/ctrlParam  ${kp_s} 0.0 0.0
fix_modify wall4_servo servo/ctrlParam  ${kp_s} 0.0 0.0

label      loop
variable   a loop 1

	# APPLING LOADING
	fix_modify	twall_servo servo/target_val -${f_max}
	fix_modify	bwall_servo servo/target_val ${f_max}
	run		${load_step} 

	# APPLING UNLOADING
	fix_modify	twall_servo servo/ctrlParam  ${kp_t1} 0.0 0.0
	fix_modify 	bwall_servo servo/ctrlParam  ${kp_t1} 0.0 0.0
	
	fix_modify	twall_servo servo/target_val -${f_min5}
	fix_modify	bwall_servo servo/target_val ${f_min5}
	run		${unload_step} 
	
	fix_modify	twall_servo servo/target_val -${f_min4}
	fix_modify	bwall_servo servo/target_val ${f_min4}
	run		${unload_step}
	
	fix_modify	twall_servo servo/target_val -${f_min3}
	fix_modify	bwall_servo servo/target_val ${f_min3}
	run		${unload_step}
	
	fix_modify	twall_servo servo/target_val -${f_min2}
	fix_modify	bwall_servo servo/target_val ${f_min2}
	run		${unload_step}
	
	fix_modify	twall_servo servo/target_val -${f_min1}
	fix_modify	bwall_servo servo/target_val ${f_min1}
	run		${unload_step}

	# RELAX SAMPLE
	run		${relax_step} 

next	   a
jump       SELF loop


#===============================================================================================================