#Particle packing by insertion and successive growing of particles

atom_style	granular
atom_modify	map array
boundary	m m m
newton		off
#echo		both

communicate	single vel yes

units		si

region		reg block -0.05 0.05 -0.05 0.05 0. 0.15 units box
create_box	1 reg

neighbor	0.002 bin
neigh_modify	delay 0


#Material properties required for new pair styles

fix 		m1 all property/global youngsModulus peratomtype 5.e6
fix 		m2 all property/global poissonsRatio peratomtype 0.45
fix 		m3 all property/global coefficientRestitution peratomtypepair 1 0.3
fix 		m4 all property/global coefficientFriction peratomtypepair 1 0.5

#New pair style
pair_style gran model hertz tangential history #Hertzian without cohesion
pair_coeff	* *

timestep	0.00001

fix xwalls1 all wall/gran model hertz tangential history primitive type 1 xplane -0.05
fix xwalls2 all wall/gran model hertz tangential history primitive type 1 xplane +0.05
fix ywalls1 all wall/gran model hertz tangential history primitive type 1 yplane -0.05
fix ywalls2 all wall/gran model hertz tangential history primitive type 1 yplane +0.05
fix zwalls1 all wall/gran model hertz tangential history primitive type 1 zplane  0.00
fix zwalls2 all wall/gran model hertz tangential history primitive type 1 zplane  0.15


#distributions for insertion
fix		pts1 all particletemplate/sphere 15485863 atom_type 1 density constant 2500 radius constant 0.005
fix		pts2 all particletemplate/sphere 15485867 atom_type 1 density constant 2500 radius constant 0.008
fix		pdd1 all particledistribution/discrete 32452843 2 pts1 0.3 pts2 0.7

#parameters for gradually growing particle diameter
variable	alphastart equal 0.25
variable	alphatarget equal 0.67
variable	growts equal 50000
variable	growevery equal 40
variable	relaxts equal 20000

#region and insertion
group		nve_group region reg

#particle insertion
fix		ins nve_group insert/pack seed 32452867 distributiontemplate pdd1 &
			maxattempt 200 insert_every once overlapcheck yes all_in yes vel constant 0. 0. 0. &
			region reg volumefraction_region ${alphastart}

#apply nve integration to all particles that are inserted as single particles
fix		integr nve_group nve/sphere

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

#insert the first particles
run		1
dump		dmp all custom/vtk 350 post/packing_*.vtk id type type x y z ix iy iz vx vy vz fx fy fz omegax omegay omegaz radius 
unfix		ins

#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
run		${growts}

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


