height =  0.40e-3
--nHeight = 8.0*2^0
dx = 5.0e-5/2^0 --height/nHeight
l_h = 3

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

--dx = height/nHeight
nHeight = math.ceil(height/dx)
nLength = nHeight*l_h 
length = nLength*dx
--length = 2.2
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 = 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-9
--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*1e4
diff_H2O_Cl = 2.097388e-9*1e4
diff_Na_Cl = 2.95407e-11*1e4

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

--flow parameters
nu_phy = 1e-6 --m^2/s
rho0_p = 1000.0
rho0_L = 1.0
cs2 = 1.0/3.0
Re = 50*1e-0
u_in_phy = Re*nu_phy/height
u_mean_phy = 2.0*u_in_phy/3.0

u_mean_L = u_mean_phy*dt/dx
u_in_L = 3.0*u_mean_L/2.0

ref_press = 1.0--( rho0_L*cs2 ) * ( rho0_p * (dx*dx) / (dt*dt) )
--press_drop = 1e-4
press_drop = 8*u_in_phy*rho0_p*nu_phy*length/height^2
press_grad = press_drop/length

function u_inflow(x,y,z) 
  return 4.0*u_in_phy*(height*0.5-y)*(height*0.5+y)/height^2.0
end

u_in_phy = u_inflow(0.0, 0.0 ,0.0)
u_mean_phy = 2.0*u_in_phy/3.0

Re = u_in_phy*height/nu_phy

nL = math.ceil(height/dx)
Re_check = u_in_phy*height/nu_phy

function omega_spatial(x,y,z)
  if x<1.0 then
    return x
  else
    return 1.0
  end
end

-- Analytical solutions
function wssA(x, y, z)
  return nu_phy * 2 * math.abs(y) * u_in_phy / (height*0.5)^2.0
end
function strainRateA(x, y, z)
  return -1*(y) * u_in_phy / R / R
end
function pressA(x,y,z)
  return ref_press+press_drop*(length-x)/length
end
