-- all units are in physical units
-- length (m)
-- mesh info
-- length of the diffusion tube
height =  1.0e-6 --m
nHeight = 100*1.00
dx = height/nHeight
l_h = 5.00

bc_type = 'wall'
bc_type = 'moments_BB'
bc_type = 'moments'
bc_type = 'BounceBack'
if bc_type == 'moments' then
  dx_offset = dx
end
if bc_type == 'wall' then
  dx_offset = dx
end
if bc_type == 'BounceBack' then
  dx_offset = 0.0
  nHeight = nHeight --+ 1  
end
if bc_type == 'moments_BB' then
  dx_offset = 0.0
  nHeight = nHeight + 0 
end

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
offset = 0--height/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 = 1.0e-8
--diff_H2O_Na = 1.31570699e-9
--diff_H2O_Cl = 2.097388e-9
--diff_Na_Cl = 2.95407e-11

-- fake diffusivities [m^2/s]
diff_diag = 1.0e-8
diff_H2O_Na = 1.31570699e-9
diff_H2O_Cl = 2.097388e-9
diff_Na_Cl = 2.95407e-11

resi_diag = 1.0/diff_diag
resi_H2O_Na = 1.0/diff_H2O_Na
-- lattice resistivity
resi_latt = 6.0
-- time step size 
dt = (resi_diag/resi_latt)*dx*dx
tmax_p = 10.0 --sec
tmax = tmax_p/dt

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

-- density
rho0 = 1.0 -- kg/m^3
-- shear viscosity
--nu = 1e-6 --m^2/s
omega_diff = 2.0
omega_kine = 2.0
nu = (1.0/(3.0*omega_kine))*dx^2/dt

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

originX = stripe_cen
originY = height/2.0
originZ = 00.0
halfwidth = width
amplitude = 0.1
background = tau

function ic_1Dgauss_pulse(x, y, z)
  return background+amplitude*math.exp(-math.log(2.)/(halfwidth^2)*( x - originX )^2)
end
function ic_2Dgauss_pulse(x, y, z)
  return background+amplitude*math.exp(-math.log(2.)/(halfwidth^2)*(( x - originX )^2+( y - originY )^2))
end


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

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

function BC_velocity(x,y,z,t)
  return deltaP * height_new*height_new*(1.0-((y-height_new/2.0)/(height_new/2.0))^2)/(length*1.0*nu)
  --return 0.0
end
u_in_phy = IC_velocity(0.0,height_new/2.0,0.0,0.0)
u_in_L = IC_velocity(0.0,height_new/2.0,0.0,0.0)*dt/dx

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

function moleflux_H2O(x,y,z,t)
  return {IC_H2O(x,y,z)*BC_velocity(x,y,z)*1000.0,0.0,0.0}
end


function moleflux_Na(x,y,z,t)
  return {IC_Na(x,y,z)*BC_velocity(x,y,z)*1000.0,0.0,0.0}
end


