-- all units are in physical units
-- length (m)
-- mesh info
-- length of the diffusion tube
height =  0.4e-3
nHeight = math.ceil(64)
dx = height/nHeight

l_h = 5
nLength = nHeight*l_h 
length = nLength*dx
--length = 7e-3
--nLength = math.ceil(length/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.10e-3
Dia = radius*2.0

-- simulation params
rho0_p = 1025 --kg/m^3
-- density
rho0_l = 1.0 
--dynamic viscosity
mu = 1.08e-3 --Pa s
--kinematic viscosity
nu_phy = mu / rho0_p  --m^2/s

-- molecular weights
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_Na = 1.0
charge_nr_Cl = -1.0

-- diffusivities
-- diagonals are not important
-- largest diffusivities
diff_diag = 1.0e-8
diff_H2O_Na =  diff_diag
diff_H2O_Cl = diff_diag

omega = 1.8--os.getenv('omega')
diff_L = (1.0/omega-0.5)/3.0
dt = diff_L*dx^2/diff_diag

-- max velocity
u_max_phy = 1e-2--m/3

tmax_p = 0.5
tmax = tmax_p/dt
t_ramp = 0.025

--concentration of sea water
--salinity of sea water = 35
--http://en.wikipedia.org/wiki/Seawater
conc_Na = 2e-2 --mol/m^3
conc_Cl = 2e-2 --mol/m^3
mole_density = conc_Na+conc_Cl
--print(conc_H2O,conc_Na,conc_Cl,moleDens0)
---- mole fraction
moleFrac_Na = conc_Na/mole_density
moleFrac_Cl = conc_Cl/mole_density
--print(moleFrac_H2O,moleFrac_Na,moleFrac_Cl)

--moleDens0 = 100  -- mol/m^3
--moleFrac_Na = 0.009
--moleFrac_Cl = 0.009
----moleFrac_Na = 1/3.
----moleFrac_Cl = 1/3.
--moleFrac_H2O = 1.0 - (moleFrac_Na+moleFrac_Cl)
--print(moleFrac_H2O,moleFrac_Na,moleFrac_Cl)
---------------------------------------
-- Initial condition for taylor dispersion
---------------------------------------
--u_max = deltaP*height*height/length*1.0*nu
--smallness parameter
tau = 0.0001
--sharpness
k = 2.0
stripe_cen = length/8.
width = stripe_cen*2.
chi = - math.log(1e-12)

function IC_velocity(x,y,z)
  return 4.0*u_max_phy*y*(height-y)/height^2.0
end

function velocity_fluid(x,y,z,t)
  return {4.0*u_max_phy*y*(height-y)/height^2.0,0.0,0.0}
end

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

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

velocity = { predefined = 'combined',
             temporal = {predefined='smooth', min_factor = 0.0, 
                          max_factor=1.0, from_time=0, to_time=t_ramp}, 
             spatial = BC_velocity
	    }

function electric_force(x,y,z,t)
  if t>t_ramp then
    return {0.0,1e-4,0.0}
  else
    return {0.0,0.0,0.0}
  end
end
