TITLE 'Pore Interconnectivity'

!---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
COORDINATES
CARTESIAN2('ro','Y')  { coordinate system, 1D,2D,3D, etc }
!---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
VARIABLES        { system variables }
theta
!---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
SELECT         { method controls }
STAGES = 50
!---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
DEFINITIONS    { parameter definitions }

!START INPUT PARAMETERS//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
	num_stages=50
	eta =0.5
	 L =  19
	 h = 0.020
	 P_in = 0.02
	P_out = 0
	K_min = 1e-5
	K_max = 1.1

!END INPUT PARAMETERS///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

x=L*ro/2 
z=h*Y 
P=theta*(P_out-P_in)+P_in

K = K_min+(stage-1)*((K_max-K_min)/(num_stages-1)) {Varying K from K_min to K_max}
k_z = 200 {Currently this parameter is not used in the program.}
k_x = K*k_z / ((2*h/L)^2) {This parameter is not currently used.}

!Used to define the boundary conditions:
blockage_start_bottom = -eta/L 
blockage_start_top = eta/L 

!should be unitless from 0 to 1.
norm_flux = SURF_INTEGRAL(dro(dY(theta)), 'F1')

V_x=-k_x*(2/L)*((P_out-P_in))*dro(theta)+P_in
V_z=-k_z*(1/h)*((P_out-P_in))*dY(theta)+P_in	

! INITIAL VALUES
!---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

EQUATIONS        { PDE's, one for each variable }

theta: dro(dro(theta)) + (1/K)*dY(dY(theta))=0 {non-dimensional version of above}

!---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
! CONSTRAINTS    { Integral constraints }
!---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
!EXTRUSION K=0.001,2
!---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
BOUNDARIES
  Region 'domain'
MESH_SPACING=.2

    Start 'outer' (-1,0)
      {walking around the sample starting at the bottom left corner which is unblocked }
    	value(theta)=1    line to (blockage_start_bottom,0)
		natural(theta)=0			line to (1,0)
		natural(theta)=0			line to (1,1)
		value(theta)=0		line to (blockage_start_top, 1)
		natural(theta)=0			line to (-1,1)
		natural(theta)=0			line to close
Feature Start "F1" (-1,0) line to (1,0)

!---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
! TIME 0 TO 1    { if time dependent }
!---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
MONITORS         { show progress }
!---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
PLOTS            { save result displays }
CONTOUR(theta) painted

HISTORY (norm_flux) VS K

END