-- all units are in physical units
-- length (m)
-- mesh info
-- length of the diffusion tube
height =  0.41
nHeight = 50*8
l_h = 5
nLength = nHeight*l_h 
dx = height/nHeight
length = nLength*dx
nLength_bnd = nLength+2
level = math.ceil(math.log(nLength_bnd)/math.log(2))
length_bnd = (2^level)*dx

dx_half = dx*0.5
zpos = dx_half
sph_pos = {0.2,0.2,zpos}
radius = 0.05
Dia = radius*2.0

-- simulation params
mH2O = 18.01528e-3 -- H20
mNa = 22.98977e-3 -- Na
mCl = 35.4527e-3  -- Cl
m_min = math.min(mH2O,mNa,mCl)
-- diffusivities
-- diagonals are not important
-- largest diffusivities
--diff_diag = 3.0e-8
--diff_H2O_Na = 1.31570699e-9
--diff_H2O_Cl = 2.097388e-9
--diff_Na_Cl = 2.95407e-11

-- fake diffusivities
diff_diag = 1.0e-2
diff_H2O_Na = 1.31570699e-6
diff_H2O_Cl = 2.097388e-3
diff_Na_Cl = 2.95407e-3

tmax_p = 50.0
-- density
rho0 = 1.0 -- kg/m^3
nL = math.ceil(Dia/dx)
omega_diff = 2.0
nu_l = (1.0/(3.0*omega_diff))

--Re = 20
---- shear viscosity
--nu = 1e-3--diff_H2O_Na
--u_mean = Re*nu/Dia
--u_max = 3.0*u_mean/2.0
--
--dt = nu_l*dx^2/nu
--u_max_l = u_max*dt/dx
--u_mean_l = 2.0*u_max_l/3.0
--nu = (1.0/(3.0*omega_diff))*dx^2/dt

dt = 1e-3
u_max_l = 0.03
u_mean_l = 2.0*u_max_l/3.0
nu = nu_l*dx^2/dt
u_max = u_max_l*dx/dt
u_mean = u_mean_l*dx/dt



--tmax = 1000
-- time step size 
--dt = 1e-2--tmax_p/tmax
tmax = tmax_p/dt

-- smallness parameters to avoid division by zero
delta = 0.0001
pressure = 99.4e3 --Pa (N/m^2)
-- Temperature in vapor phase
temp = 328.5 --K

---------------------------------------
-- Boundary condition for stefan tube
---------------------------------------
-- At Inlet, x=0:
moleFrac_Ace_In = 0.319
moleFrac_Met_In = 0.528
moleFrac_Air_In = 1.0 - ( moleFrac_Ace_In + moleFrac_Met_In )

-- At outlet, x=length:
moleFrac_Ace_Out = delta
moleFrac_Met_Out = delta
moleFrac_Air_Out = 1.0 - ( moleFrac_Ace_Out + moleFrac_Met_Out )

---------------------------------------
-- Initial condition for taylor dispersion
---------------------------------------
--c_h = height
--deltaP = 1e-2
--umax = deltaP * height*height/(length*1.0*nu)
--smallness parameter
tau = 0.0001
--sharpness
k = 2.0
stripe_cen = length/16.
width = stripe_cen*1.
chi = - math.log(1e-12)

function IC_pressure(x,y,z)
  return rho0*(-deltaP*x/length + deltaP)
end

function IC_velocity(x,y,z)
  --return deltaP * height*height*(1.0-((y-height/2.0)/(height/2.0))^2)/(length*1.0*nu)
  return 4.0*u_max*y*(height-y)/height^2.0
end

function BC_velocity(x,y,z,t)
  return IC_velocity(x,y,z)
end

function IC_Na(x,y,z)
  return tau + 0.1 * math.exp(-chi*((x-stripe_cen)/(width/2.0))^(2.0*k))
end

function BC_Na(x,y,z,t)
  return IC_Na(x,y,z)
end

function IC_H2O(x,y,z)
  na = IC_Na(x,y,z)
  return 1.0 - na
end

function BC_H2O(x,y,z,t)
  na = IC_Na(x,y,z)
  return 1.0- na
end



-- boundary condition
bc = {   
--   { kind = 'wall',
--     label = 'east'},
--   { kind = 'wall',
--     label = 'west'},
}
