--mesh parameters
--height of the spacer channel
h_ch = 0.4e-3 --m
--diameter of the spacer filament
d_f = h_ch*0.5
--radius of the filament
r_f = d_f*0.5
--mesh_length to height ratio
lm_hsp = 2.5
--mesh length or distance between spacer filament
l_m = h_ch*lm_hsp
--number of spacer filament in x-direction
nFilament = 2.0
--length of the channel
l_ch = l_m*nFilament
--width of the channel
w_ch = l_m*2.0
--Height of the channel in lattice unit/number of elements along the height
nHeight = 20
--discretization size
dx = h_ch/nHeight
dx_half = dx*0.5
--Find refinement level required for this dx
nLength = math.ceil(l_ch/dx)
nWidth = math.ceil(w_ch/dx) 
nLength_bnd = nLength+2
level = math.ceil(math.log(nLength_bnd)/math.log(2))
length_bnd = (2^level)*dx

-- Flow parameter
--compute hydralic diameter
voidage = 1.0 - (math.pi*d_f^2.0/(2.0*l_m*h_ch))
d_h = 4.0*voidage/(2./h_ch + (1.0-voidage)*(4.0/d_f))
--density of pure water
--taken from "Perry Handbook of chemical engineering"
rho0_p = 998.2071 --kg/m^3
--dynamic viscosity
mu = 1.08e-3 --Pa s
--kinematic viscosity
nu_phy = mu / rho0_p  --m^2/s

--nElements on filement diameter
nL = math.ceil(d_f/dx)
--nElements on hydralic diameter
nL_h = math.ceil(d_h/dx)

u_mean_phy = 1e-3 --m/s
u_max_phy = u_mean_phy*3.0/2.0 --m/s
Re = u_mean_phy*d_h/nu_phy

--omega
omega_diff = 2.0
omega_kine = 0.1
nu_L = (1.0/(3.0*omega_kine))
dt = nu_L*dx^2./nu_phy
u_max_L = u_max_phy*dt/dx
u_mean_L = 2.0*u_max_L/3.0

tmax_p = 1.0
tmax = tmax_p/dt

--u_max_L = 0.05
--u_mean_L = 2.0*u_max_L/3.0

--dt = u_max_L*dx/u_max_phy
--nu_L = nu_phy*dt/dx^2
Re_l = u_mean_L*nL_h/nu_L
----Reynolds number
--Re = d_h*u_mean_phy/nu_phy
--Re_l = nL_h*u_mean_L/nu_L

seedPoint = {l_m,d_f,w_ch/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)

-- specific charge
charge_nr_H2O = 0.0
charge_nr_Na = 1.0
charge_nr_Cl = -1.0

-- 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-9
diff_H2O_Na = 1.31570699e-9
diff_H2O_Cl = 2.097388e-9
diff_Na_Cl = 2.95407e-11

-- 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

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

function IC_velocity(x,y,z)
  return 4.0*u_max_phy*y*(h_ch-y)/h_ch^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_Cl(x,y,z)
  --return tau + 0.1 * math.exp(-chi*((x-stripe_cen)/(width/2.0))^(2.0*k))
  return IC_Na(x,y,z)
end

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


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

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




