-- all units are in physical units
-- length (m)
-- mesh info
-- length of the diffusion tube
length = 0.238 --m
nLength = 60*1
nWidth = 2
nDepth = 2
dx = length/nLength
nLength_bnd = nLength + 2
level = math.ceil(math.log(nLength_bnd)/math.log(2))
length_bnd = (2^level)*dx

-- simulation params
tmax_p = 100000.00
--tmax = 1000
-- time step size 
dt = 1e+1--tmax_p/tmax
tmax = tmax_p/dt

-- Stefan tube testcase paraClers
-- REF: Multicomponent Mass Transfer - R. Taylor, R. Krishna. Pg. 21
-- molecular weights(kg/mol)
mH2O = 18.01528e-3 -- H20
mNa = 22.98977e-3 -- Na
mCl = 35.4527e-3  -- Cl
-- binary diffusivities in m^2/s
diff_diag = 1.0e-2
diff_H2O_Na = 1.31570699e-9
diff_H2O_Cl = 2.097388e-9
diff_Na_Cl = 2.95407e-13

-- pressure in vapor phase
pressure = 99.4e3 --Pa (N/m^2)
-- Temperature in vapor phase
temp = 328.5 --K

-- density
rho0 = 1000.0 -- kg/m^3

-- smallness paraClers to avoid division by zero
delta = 0.0001

---------------------------------------
-- Boundary condition for stefan tube
---------------------------------------
-- At Inlet, x=0:
moleFrac_Na_In = 0.01
moleFrac_Cl_In = 0.01
moleFrac_H2O_In = 1.0 - ( moleFrac_Na_In + moleFrac_Cl_In )

-- At outlet, x=length:
moleFrac_Na_Out = delta
moleFrac_Cl_Out = delta
moleFrac_H2O_Out = 1.0 - ( moleFrac_Na_Out + moleFrac_Cl_Out )

---------------------------------------
-- Initial condition for stefan tube
---------------------------------------
function IC_Na_Linear(x,y,z)
  return moleFrac_Na_In * (1.0 - x/length) + moleFrac_Na_Out
end

function IC_Na_Smooth(x,y,z)
  return moleFrac_Na_In * 0.5 * (1.0 - math.tanh((x-length/2.0)/dx)) + delta
end

function IC_Cl_Linear(x,y,z)
  return moleFrac_Cl_In * (1.0 - x/length) + moleFrac_Cl_Out
end

function IC_Cl_Smooth(x,y,z)
  return moleFrac_Cl_In * 0.5 * (1.0 - math.tanh((x-length/2.0)/dx)) + delta
end

function IC_H2O_Linear(x,y,z)
  ic_Na = IC_Na_Linear(x,y,z)
  ic_Cl = IC_Cl_Linear(x,y,z) 
  return 1.0 - ( ic_Na + ic_Cl )
end

function IC_H2O_Smooth(x,y,z)
  ic_Na = IC_Na_Smooth(x,y,z) 
  ic_Cl = IC_Cl_Smooth(x,y,z) 
  return 1.0 - ( ic_Na + ic_Cl )
end


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