#Particle packing by insertion by gravity 
# DST packing, with respect to Ni et al. 2000, Leighton Buzzard sand, Shearbox 60 mm x 60 mm x 20 mm 
# Constant settings
atom_style		granular
atom_modify		map array
boundary		f f f  # will change to f f f boundary to save time 
newton			off
#echo			both
communicate		single vel yes
units			si

region			reg block -0.03 0.03 -0.02 0.02 0.0 0.015 units box 
create_box		2 reg

neighbor		0.00064 bin
neigh_modify		delay 0


#Material properties required for new pair styles

fix 			m1 all property/global youngsModulus peratomtype 2.6e8 2.6e8 
fix 			m2 all property/global poissonsRatio peratomtype 0.3 0.3
fix 			m3 all property/global coefficientRestitution peratomtypepair 2 0.2 0.2 0.2 0.2
fix 			m4 all property/global coefficientFriction peratomtypepair 2 0.0 0.0 0.0 0.0
fix 			m5 all property/global coefficientRollingFriction peratomtypepair 2 0.0 0.0 0.0 0.0
#fix 			m6 all property/global characteristicVelocity scalar 2
fix 			m7 all property/global cohesionEnergyDensity peratomtypepair 2 0 0 0 0

#New pair style
pair_style 		gran model hertz tangential history cohesion sjkr rolling_friction epsd2 
pair_coeff		* *

timestep		1e-6 # 17% tcrit

## Primitive walls to minimize the use of meshed walls, to save time 
fix  			xwalls1 all wall/gran model hertz tangential history primitive type 2 xplane -0.015
fix 			xwalls2 all wall/gran model hertz tangential history primitive type 2 xplane  0.015
fix 			ywalls1 all wall/gran model hertz tangential history primitive type 2 yplane -0.015
fix 			ywalls2 all wall/gran model hertz tangential history primitive type 2 yplane  0.015
fix 			zwalls1 all wall/gran model hertz tangential history primitive type 2 zplane  0.00
fix 			zwalls3 all wall/gran model hertz tangential history primitive type 2 zplane  0.01


#distributions for insertion
# 15485863, 15485867, 32452843, 32452867, 49979687, 49979693, 67867967, 67867979, 86028121, 86028157, 98653427, 96835421
fix			pts1 all particletemplate/sphere 15485863 atom_type 1 density constant 2650 radius constant 0.00064
fix			pts2 all particletemplate/sphere 15485867 atom_type 1 density constant 2650 radius constant 0.00078
fix			pts3 all particletemplate/sphere 32452843 atom_type 1 density constant 2650 radius constant 0.00081
fix			pts4 all particletemplate/sphere 32452867 atom_type 1 density constant 2650 radius constant 0.001
fix			pts5 all particletemplate/sphere 49979687 atom_type 1 density constant 2650 radius constant 0.0012

fix			pdd1 all particledistribution/discrete 49979693 5 pts1 0.1 pts2 0.4 pts3 0.1 pts4 0.2 pts5 0.2 

#region and insertion
region			reg_ins block -0.015 0.015 -0.015 0.015 0.0 0.01 units box

#parameters for gradually growing particle diameter
variable		alphastart equal 0.4
variable		alphatarget equal 0.7
variable		growts equal 100000
variable		growevery equal 40
variable		relaxts equal 100000
#particle insertion
			
fix 			ins all insert/pack seed 67867967 distributiontemplate pdd1 maxattempt 1000 insert_every once overlapcheck yes all_in no vel constant 0. 0. -1 &
			region reg_ins volumefraction_region ${alphastart}


#apply nve integration to all particles that are inserted as single particles
fix			integr all nve/sphere
fix			gravi all gravity 9.81 vector 0.0 0.0 -1.0


#output settings, include total thermal energy
compute		1 all erotate/sphere
thermo_style		custom step atoms ke c_1 vol
thermo			10000
thermo_modify		lost ignore norm no

#insert the first particles
fix			ts all check/timestep/gran 1 0.2 0.2			# for loop, fix/time/check should be placed after fix/couple/cfd
run			1
unfix			ts

dump			dmp all custom/vtk 10000 post1/Particle/packing_*.vtk id type type x y z ix iy iz vx vy vz fx fy fz omegax omegay omegaz radius 

# Run insertion 
run			100000
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     		100000 						# Comment it while reading restart file
#unfix			grow

#calculate grow rate
variable		Rgrowrate equal (${alphatarget}/${alphastart})^(${growevery}/(3.*${growts}))
print			"The radius grow rate is ${Rgrowrate}"

#do the diameter grow
compute 		rad all property/atom radius

variable		dgrown atom ${Rgrowrate}*2.*c_rad
fix			grow all adapt ${growevery} atom diameter v_dgrown
run			${growts}

#let the packing relax
unfix	      	  	grow
run			${relaxts}

unfix			xwalls1
unfix			xwalls2
unfix			ywalls1
unfix			ywalls2
unfix			zwalls1
unfix			zwalls3 
write_restart 		liggghts.packing_gravity1