# Khane 2016 Simulation

# RADIUS : 0.00635 m (0.25 in)
# DIAMETER : 0.0127 m (0.50 in)

# ========================= Initialization ========================== #
# =================================================================== #

atom_style 	sphere
atom_modify 	map array
units 		si
newton 		off
dimension 	3
boundary 	f f f
processors 	2 2 4
communicate 	single vel yes

# Define pebble dimensions
variable 	pdiam equal 0.0127
variable 	prad equal 0.00635
variable 	diameter_margin equal 1 # +/- this many diameters on insertion zones
variable 	insert_margin equal ${pdiam}*${diameter_margin}

neighbor 	${prad} bin
neigh_modify	delay 0

# ================== Atom Definition and Insertion ================== #
# =================================================================== #

region 		box_domain block -0.16 0.16 -0.16 0.16 -0.1 0.5 units box # Try changing later for cylinder
create_box 	2 box_domain

fix  		pts1 all particletemplate/sphere 15485863 atom_type 1 density constant 1540 radius constant ${prad}
fix  		pdd1 all particledistribution/discrete 49979687 1 pts1 1.0

# Tracer group
group		tracers type 2

fix		pts2 tracers particletemplate/sphere 15485863 atom_type 2 density constant 1540 radius constant ${prad}
fix 		pdd2 tracers particledistribution/discrete 49979687 1 pts2 1.0

# Define insertion region, include cushioning (~1 diam) to prevent intersection with walls
		 ##    ID      style  dim c1 c2 radius  lo   hi
region 		insert_zone cylinder  z   0  0  0.135   0.331  0.38 units box
region 		small_insert_zone cylinder  z   0  0  0.032   0.342  0.38 units box

# Define four radial insertion regions (r/R = 0, 0.33, 0.67, 0.92), include cushioning (~1 diam) to prevent intersection with walls
		 ##    ID      style  dim c1 c2 radius  lo   hi
	# r/R = 0
region 		rR_0 cylinder  z   0  0  ${insert_margin}   0.331  0.38 units box

	# r/R = 0.33
region 		rR_33_inner cylinder  z   0  0  0.0495-${insert_margin}   0.331  0.38 units box side out
#region 		rR_33_inner cylinder  z   0  0  0.0368   0.331  0.38 units box side out
region 		rR_33_outer cylinder  z   0  0  0.0495+${insert_margin}   0.331  0.38 units box side in
region 		rR_33 intersect 2 rR_33_inner rR_33_outer

	# r/R = 0.67
region 		rR_67_inner cylinder  z   0  0  0.1005-${insert_margin}   0.331  0.38 units box side out
region 		rR_67_outer cylinder  z   0  0  0.1005+${insert_margin}   0.331  0.38 units box side in
region 		rR_67 intersect 2 rR_67_inner rR_67_outer

	# r/R = 0.92
region 		rR_92_inner cylinder  z   0  0  0.1373-${insert_margin}  0.331  0.38 units box side out
region 		rR_92_outer cylinder  z   0  0  0.1373+${insert_margin}  0.331  0.38 units box side in
region 		rR_92 intersect 2 rR_92_inner rR_92_outer

# Insert initial load of pebbles
# Full = 12,000, Partial = 10,000 (leaves insertion region open)
fix 		ins1 all insert/rate/region seed 32452843 distributiontemplate pdd1 nparticles 11000 particlerate 20000 insert_every 10000 overlapcheck yes vel constant 0. 0. -2.0 region insert_zone

# ============================ Settings ============================= #
# =================================================================== #

# Contact model
pair_style 	gran model hooke/stiffness tangential history tangential_damping on absolute_damping on limitForce off
pair_coeff 	* *

# Numerical scheme
#fix 		integrate all nve/sphere
#timestep 	0.000001
timestep 	0.000002
hard_particles	yes # glass has large YM

# Glass properties
variable 	CF equal 0.2353
variable	w1 atom 0.0
variable 	w2 atom 0.0
variable 	YM equal 2.4e10
variable 	PR equal 0.25
variable 	kn_val equal 76200000
variable 	kt_val equal 30480000
variable 	gamman_val equal 69.237
variable 	gammat_val equal 43.789

# Number of tracers
variable	N_tracer equal 600/4

# Material properties for hooke/stiffness with absolute damping on
fix  		m1 all property/global coefficientFriction peratomtypepair 2 ${CF} ${CF} ${CF} ${CF}
fix 		m2 all property/global youngsModulus peratomtype ${YM} ${YM}
fix 		m3 all property/global poissonsRatio peratomtype ${PR} ${PR}
fix		p1 all property/global kn peratomtypepair 2 ${kn_val} ${kn_val} ${kn_val} ${kn_val}																										
fix		p2 all property/global kt peratomtypepair 2 ${kt_val} ${kt_val} ${kt_val} ${kt_val}																										
fix		p3 all property/global gamman_abs peratomtypepair 2 ${gamman_val} ${gamman_val} ${gamman_val} ${gamman_val}																									
fix		p4 all property/global gammat_abs peratomtypepair 2 ${gammat_val} ${gammat_val} ${gammat_val} ${gammat_val}

# Use the following to check if timestep is adequate based on Young's Modulus
#fix 		ts_check all check/timestep/gran 1000 0.1 0.1	 

# Gravity
#fix 		grav all gravity 9.81 vector 0.0 0.0 -1.0

# Import mesh from cad (make sure to rotate on x-axis by 90deg):
fix 		pbr all mesh/surface file mesh/Khane_Core_Geometry_OUTLET.stl type 1 heal auto_remove_duplicates scale 0.001 move -0.15 0 -0.15 rotate axis 1 0 0 angle 90 curvature_tolerant yes 

# Use the imported mesh as granular wall
fix 		geometry all wall/gran model hooke/stiffness tangential history mesh n_meshes 1 meshes pbr tangential_damping on absolute_damping on limitForce off

# Let the pebbles fill up against stopper before flowing out
fix 		stopper all wall/gran model hooke/stiffness tangential history primitive type 1 zplane 0.0 tangential_damping on absolute_damping on limitForce off

# ========================== Run and Dump =========================== #
# =================================================================== #

# Initialize
#run    		1

# Produce figures for vessel filling
#dump   		dmp1 all custom/vtk 10000 post/filling/khane_filling_*.vtk id type type x y z ix iy iz vx vy vz radius
#dump 		d1 all image 10000 post/filling/khane_filling_*.jpg type type adiam ${pdiam} view 60 300 zoom 1.2 box yes 0.03 axes yes 0.35 0.03

# Fill vessel (run a little extra long to settle)
#run 		550000

unfix 		ins1
#undump 		dmp1
#undump		d1

# Add tracer pebbles to top (r/R = 0, 0.33, 0.67, 0.92)
fix 		ins2_0 tracers insert/rate/region seed 11657 distributiontemplate pdd2 nparticles ${N_tracer} particlerate 7500 insert_every 1000 overlapcheck yes all_in no vel constant 0. 0. 0. region rR_0
fix 		ins2_33 tracers insert/rate/region seed 11657 distributiontemplate pdd2 nparticles ${N_tracer} particlerate 7500 insert_every 1000 overlapcheck yes all_in no vel constant 0. 0. 0. region rR_33
fix 		ins2_67 tracers insert/rate/region seed 11657 distributiontemplate pdd2 nparticles ${N_tracer} particlerate 7500 insert_every 1000 overlapcheck yes all_in no vel constant 0. 0. 0. region rR_67
fix 		ins2_92 tracers insert/rate/region seed 11657 distributiontemplate pdd2 nparticles ${N_tracer} particlerate 7500 insert_every 1000 overlapcheck yes all_in no vel constant 0. 0. 0. region rR_92

run 		10000

unfix 		ins2_0
unfix 		ins2_33
unfix 		ins2_67
unfix 		ins2_92

# Let bed settle
#run 		300000

# Remove bottom to simulate outflow
unfix 		stopper

# Add pebbles as we go (INF gives us as many as we may need)
#fix 		ins3 all insert/rate/region seed 11657 distributiontemplate pdd1 nparticles INF particlerate 20000 insert_every 10000 overlapcheck yes vel constant 0. 0. -2.0 region small_insert_zone

# Calculate residence times (define region + define tracer variable)
		 ##    ID      style  dim c1 c2 radius  lo   hi
region 		main_cylinder cylinder z 0 0 0.15 0.08 0.38
fix 		rttracker all property/atom/regiontracer/time add_region main_cylinder # rttracker = [ RT in entire simulation , RT in main_cylinder ]
#fix 		rttracker all property/atom/regiontracer/time

# Produce figures/data files for quasi-steady flow
dump   		dmp2 all custom/vtk 15000 post/khane_bulk_*.vtk id type type x y z vx vy vz radius f_rttracker[1] f_rttracker[2]
dump   		dmp3 tracers custom/vtk 15000 post/khane_tracers_*.vtk id type type x y z vx vy vz radius f_rttracker[1] f_rttracker[2]
dump 		d2 tracers image 15000 post/khane_tracers_*.jpg type type adiam ${pdiam} view 60 300 zoom 1.2 axes yes 0.35 0.03
dump		dat1 tracers custom 10000 post/khane_test.dat id x y z radius f_rttracker[1] f_rttracker[2] v_w1 v_w2
#dump 		csv1 tracers custom 15000 post/rt_data.csv id f_rttracer f_rttracker[1] f_rttracker[2]
#dump 		csv2 all custom 1000 post/rt_data.csv id f_rttracer f_rttracker[1] f_rttracker[2] x y z

# Begin full simulation
#run		100000
#run 		5000000
#run 		6300000
#run 		7000000
#run 		10000000

##### TO-DO
# PEBBLE INSERTION
# FILL - INSERT TRACERS - DRAIN&FILL STEADY STATE