###################################Header for General commands#######################################################

#particle packing by insertion and successive growing of particles

atom_style		granular          #Simulation of particles
atom_modify            map array
boundary               p f f             #fixed boundaries -> particles will be deleted if leaving the simulation box and periodic boundaries in x axis
newton                 off               #default
communicate            single vel yes
units                  si





###################################System variables##################################################################

#Definition of boundaries
variable xmin_reg equal -1e-5		
variable xmax_reg equal 3.0001e-3		

variable ymin_reg equal -0.05001e-3
variable ymax_reg equal 8e-4

variable zmin_reg equal -1e-8
variable zmax_reg equal 3.0001e-3




variable xmin_inbox equal 0		
variable xmax_inbox equal 3e-3		

variable ymin_inbox equal 0
variable ymax_inbox equal 7.5e-4

variable zmin_inbox equal 0
variable zmax_inbox equal 3e-3 




variable xmin_leftwall equal 1e-8		
variable xmax_leftwall equal 2.999e-3		

variable ymin_leftwall equal 7.499e-4
variable ymax_leftwall equal 7.99e-4

variable zmin_leftwall equal 1e-8
variable zmax_leftwall equal 2.99e-3



variable xmin_rightwall equal 1e-8		
variable xmax_rightwall equal 2.999e-3		

variable ymin_rightwall equal -0.049e-3
variable ymax_rightwall equal -1e-8

variable zmin_rightwall equal 1e-8
variable zmax_rightwall equal 2.999e-3



###################################Definition of simulationbox#######################################################

region		         reg block ${xmin_reg} ${xmax_reg} ${ymin_reg} ${ymax_reg} ${zmin_reg} ${zmax_reg} units box
region                  inbox block ${xmin_inbox} ${xmax_inbox} ${ymin_inbox} ${ymax_inbox} ${zmin_inbox} ${zmax_inbox} units box
region                  leftwall block ${xmin_leftwall} ${xmax_leftwall} ${ymin_leftwall} ${ymax_leftwall} ${zmin_leftwall} ${zmax_leftwall} units box
region                  rightwall block ${xmin_rightwall} ${xmax_rightwall} ${ymin_rightwall} ${ymax_rightwall} ${zmin_rightwall} ${zmax_rightwall} units box
create_box              1 reg







###################################neighbor list for element interaction######################################################
neighbor                0.0005 bin
neigh_modify            delay 0             #default








###################################Definition of Material properties#################################################

#Material properties required for new pair styles

fix                      m1 all property/global youngsModulus peratomtype 7.e10
fix                      m2 all property/global poissonsRatio peratomtype 0.20 
fix                      m3 all property/global coefficientRestitution peratomtypepair 1 0.3
fix                      m4 all property/global coefficientFriction peratomtypepair 1 0.5
fix                      m5 all property/global coefficientRollingFriction peratomtypepair 1 0










#################################Definition of the contact models####################################################

#New pair style

pair_style              gran model hertz tangential history rolling_friction cdt #Hertzian without cohesion contact model 
pair_coeff              * *                    #default





###################################Generation and Loading of the Geometry .stl#######################################

fix                      ywall1 all wall/gran model hertz tangential history primitive type 1 yplane -0.05001e-3
fix                      ywall2 all wall/gran model hertz tangential history primitive type 1 yplane 8e-4
fix                      barrier1 all wall/gran model hertz tangential history primitive type 1 yplane 0 
fix                      barrier2 all wall/gran model hertz tangential history primitive type 1 yplane 7.5e-4 
fix                      xwall1 all wall/gran model hertz tangential history primitive type 1 xplane -1e-5
fix                      xwall2 all wall/gran model hertz tangential history primitive type 1 xplane 3.0001e-3 
fix                      zwall all wall/gran model hertz tangential history primitive type 1 zplane -1e-8
fix                      zwall_out all wall/gran model hertz tangential history primitive type 1 zplane 3.0001e-3 

#Definition of the timestep

timestep                 1e-9
fix                      ts_check all check/timestep/gran 1000 0.2 0.2






###################################Generation and Insertion of the particles#########################################

#distributions for insertion

fix                      pts1 all particletemplate/sphere 15485863 atom_type 1 density constant 2650 radius constant 1e-4
fix                      pts2 all particletemplate/sphere 15485867 atom_type 1 density constant 2650 radius constant 9e-5
fix                      pdd1 all particledistribution/discrete 32452843 2 pts1 1.0 pts2 0 
fix                      pdd2 all particledistribution/discrete 32452867 2 pts1 0 pts2 1.0




###################################Specific variables for current simulation########################################

#parameters for gradually growing particle diameter

variable      alphastart     equal     0.40
variable      alphatarget    equal     0.65
variable      growts         equal     350000
variable      growevery      equal     40
variable      relaxts        equal     250000


#region and insertion

group                  nve_group region reg
fix                    ins nve_group insert/stream seed 49979687 distributiontemplate pdd1 vol ${alphastart} 200 region inbox


#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                 100000
thermo_modify          lost ignore norm no
compute_modify         thermo_temp dynamic yes
print                  "Inserting Proppant Pack"


#insert the first particles

run                   1


#dump                 dmp_paraview all custom 4000 ./post/dump_forces_limited_var_7.5e-4_1e03_*.liggghts id type type x y z ix iy iz vx vy vz fx fy fz omegax omegay omegaz radius
unfix                 ins
#undump               dumpstl


#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     d_grown atom ${Rgrowrate}*2.*c_rad
fix                  grow nve_group adapt ${growevery} atom diameter d_grown
neigh_modify         every ${growevery} check no
neighbor             0.000005 bin


#run

run                  ${growts}


#let the packing relax

unfix                grow
neigh_modify         check yes
run                  ${relaxts}
variable     pack     equal      ${growts}+${relaxts}
print                "Creating Fracture Walls"
neighbor             0.0005 bin
lattice              sc 5e-5
create_atoms         1 region leftwall units box
group                leftwall_group region leftwall
set                  group leftwall_group diameter 5e-5 density 2650
create_atoms         1 region rightwall units box
group                rightwall_group region rightwall
set                  group rightwall_group diameter 5e-5 density 2650
velocity             rightwall_group zero linear
fix 2                rightwall_group nve/noforce
fix 3                rightwall_group freeze

#Granular Wall

fix                  1 leftwall_group rigid group 1 leftwall_group force * on on off torque * off off off
neigh_modify         exclude group leftwall_group leftwall_group
neigh_modify         exclude group rightwall_group rightwall_group
neigh_modify         exclude group leftwall_group rightwall_group
thermo               1000
#run 100000

unfix                barrier

#save data

variable      massleft           equal      mass(nve_group,inbox)
variable      totalmass          equal      mass(nve_group,inbox)
variable      averagevelocity    equal      vcm(nve_group,z,reg)
variable      fracturewidth      equal      xcm(leftwall_group,y,reg)-xcm(rightwall_group,y,reg)
variable      pressureleft       equal      fcm(leftwall_group,y,reg)
variable      pressureright      equal      fcm(rightwall_group,y,reg)
fix                 screen all print 100000 "${massleft}, ${totalmass}, ${averagevelocity}, ${fracturewidth}, ${pressureleft}, ${pressureright}"
fix                 extra all print 4000 "${massleft}, ${totalmass}, ${averagevelocity}, ${fracturewidth}, ${pressureleft}, ${pressureright}" file flowingMAXs_203.txt

#relax

run 100000


#compact

print              "Compacting Proppant Pack"
fix                 shmin leftwall_group addforce 0.0 0 0.0
run                 400000
fix                 damp nve_group viscous 0.189
run                 200000
unfix               damp

unfix               zwall_out
print              "Applying Fracture Flow"
fix                 gr nve_group addforce 0. 0. 8.6167e-3
run                 1000000

label               loop
print              "Entering mass loop"
variable     masspres     equal         ${massleft}
run                 40000
variable     deltam       equal         v_masspres-v_massleft
if                  $deltam=<0.00000001 then "jump in.var_width_MAXp7.5e408.6167e3 continuing"
variable     masspres     delete
variable     deltam       delete
jump               in.var_width_MAXp7.5e408.6167e3 loop
label              continuing
run                100000
print              "Done"