dimension	3
atom_style	hybrid bond meso
units		si
newton	 	on
boundary	p p p

#read_data	intestine.initial
read_restart	intestine.0

reset_timestep 0

variable	R	equal 0.025
variable	L	equal 0.6283185307180
variable	dL	equal 0.006283185307180
variable	h	equal 2.*${dL}
variable	cut	equal 1.*${dL}
variable	skin 	equal 0.3*${h} 		
variable	c	equal 0.005 # min 0.005 max 0.5
variable	rho	equal 1000
variable	alpha	equal 1. # 1. for default viscosity 
variable	nu	equal ${alpha}*${h}*${c}/10. 
variable	kB	equal 0.5
variable	f0	equal -1.e-3
print 		'kinematic viscosity is ${nu}'
variable	Nsteps		equal 50
variable	Dt		equal 0.002
variable	Nout		equal ${Nsteps}
variable	N		equal 10.
variable	Dz		equal ${L}/${N}
variable	index   	equal 0. #integer between 0 and N-1  
variable	z0		equal (v_index)*${Dz}%${L}
variable	z1		equal (v_index+1.)*${Dz}%${L}
variable	z2		equal (v_index+2.)*${Dz}%${L}
variable	vel_ind		equal 0.1 # % of L in Nsteps,  min 0.01 max 0.1

group	wall  type 1
group	fluid type 2

pair_style hybrid/overlay sph/rhosum 1 sph/taitwater soft ${cut} zero ${cut} 
pair_coeff 1 1 zero
pair_coeff 1 2 soft 2.e-6 
pair_coeff 2 2 soft 1.e-9 #artificial pressure 
pair_coeff * 2 sph/taitwater ${rho} ${c} ${alpha} ${h}
pair_coeff 2 2 sph/rhosum ${h}
bond_style	harmonic
bond_coeff	1 ${kB} ${dL}
neighbor	 ${skin} bin

variable PP atom ((v_z0<v_z1)|^((z>=v_z0)|^(z<=v_z1)))*gmask(wall)
#section that is squeezed (push) z0<z<z1
variable SS atom ((v_z1<v_z2)|^((z>v_z1)|^(z<=v_z2)))*gmask(wall)
#section that senses (sense) z1<z<z2
variable RR atom ((v_z0>v_z2)|^((z>v_z0)|^(z<v_z2)))*gmask(wall)
#the remaining section (rigid) 

group 	push    dynamic wall var PP every 1
group 	sense   dynamic wall var SS  every 1

variable rDist   atom sqrt(x^2+y^2)
variable cosQ    atom x/v_rDist
variable sinQ    atom y/v_rDist 
variable squeeze_factor  atom (v_rDist/${R})^2. 

variable FP atom ${f0}*gmask(push)  
variable Fx atom v_FP*v_cosQ*v_squeeze_factor
variable Fy atom v_FP*v_sinQ*v_squeeze_factor

###the following only because spring/self cannot be dynamic####
fix x0 wall store/state 0 x
fix y0 wall store/state 0 y
fix z0 wall store/state 0 z
fix zu wall store/state 0 zu
variable varX atom (f_x0-x)*gmask(wall)
variable varY atom (f_y0-y)*gmask(wall)
#variable varZ atom (f_z0-f_zu)*gmask(wall)
#variable varP atom sqrt(varX^2+varY^2+varZ^2)
variable KR  equal 0.1 #rigid section
variable FRx atom v_RR*${KR}*v_varX
variable FRy atom v_RR*${KR}*v_varY
#variable FRz atom v_RR*${KR}*v_varZ
#################################################################
fix 10  all     meso
fix 20	wall    viscous 0.01
fix 30  wall    addforce v_Fx v_Fy 0.
fix 40  wall    addforce v_FRx v_FRy 0.
fix 50  wall    setforce NULL NULL 0.

compute  rho_peratom all    meso/rho/atom
compute  stretch     wall   pe/atom bond
compute  reward      fluid  com
variable r_stretch   atom   gmask(sense)*((v_rDist-${R})/${R})*(v_rDist>${R}) 

thermo	10
thermo_style	 custom step time c_reward[3] v_index
thermo_modify	 norm no
dump		 my_dump all custom ${Nout} output.lammpstrj id type xs ys zs vx vy vz v_r_stretch c_rho_peratom c_stretch
dump_modify	 my_dump format float %.3e first yes

timestep	${Dt} 

variable   index   equal 1. 
run		 ${Nsteps}
reset_timestep 0

variable index_final equal ${index}
label loop
variable i loop 2
variable   Di           equal  (${vel_ind}/${Nsteps})*step
variable   index_init   equal  ${index_final}
variable   index        equal  (${index_init}+v_Di)
run		 ${Nsteps}
variable   index_final  equal  ${index_final}+${vel_ind}
reset_timestep 0
next i
jump SELF loop
label break