--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 = 64
--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 salt water
rho0_p = 1025 --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)

-- 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
fac = 1e3
diff_diag = 1.0e-9*fac
diff_H2O_Na = 1.31570699e-9*fac
diff_H2O_Cl = 2.097388e-9*fac
diff_Na_Cl = 2.95407e-11*fac

resi_diag = 1.0/diff_diag
resi_H2O_Na = 1.0/diff_H2O_Na
resi_Na_Cl = 1.0/diff_Na_Cl
resi_ref = resi_Na_Cl
resi_ref = resi_H2O_Na
-- lattice resistivity
resi_latt = 5.0
-- time step size 
dt = (resi_ref/resi_latt)*dx*dx
--print('dt =', dt)

u_max_phy = 5.0e-2 --m/3
u_mean_phy = u_max_phy*2.0/3.0 --m/s

Re = u_mean_phy*d_h/nu_phy

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

tmax_p = 0.0025
tmax = tmax_p/dt
t_ramp = 0.001

----Reynolds number
Re_l = u_mean_L*nL_h/nu_L

-- smallness parameters to avoid division by zero
pressure = 99.4e3 --Pa (N/m^2)

-------------------------------------------
-- Initial condition and boundary condition
-------------------------------------------
--concentration of sea water
--salinity of sea water = 35
--http://en.wikipedia.org/wiki/Seawater
conc_H2O = 53.6*rho0_p --mol/m^3
conc_Na = 0.469*rho0_p --mol/m^3
conc_Cl = 0.546*rho0_p --mol/m^3
moleDens0 = conc_H2O+conc_Na+conc_Cl
--print(conc_H2O,conc_Na,conc_Cl,moleDens0)
---- mole fraction
moleFrac_H2O = conc_H2O/moleDens0
moleFrac_Na = conc_Na/moleDens0
moleFrac_Cl = conc_Cl/moleDens0
--print ('molefrac', 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)

--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
  return 0.0
end

function BC_velocity(x,y,z,t)
  return {4.0*u_max_phy*y*(h_ch-y)/h_ch^2.0,0.0,0.0}
  --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))
  return tau
end

function BC_Na(x,y,z,t)
  if t>t_ramp then
    return moleFrac_Na
  else
    return tau 
  end  
end

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

function BC_Cl(x,y,z,t)
  --return IC_Na(x,y,z)
  if t>t_ramp then
    return moleFrac_Cl
  else
    return tau  
  end  
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 = BC_Na(x,y,z,t)
  cl = BC_Cl(x,y,z,t)
  return 1.0 - na - cl
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
             --spatial = {u_max_phy,0.0,0.0}
	    }

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