useObstacle = true 
height =  50.0e-9
nHeight = 256
l_h = 1
nLength = 1--nHeight*l_h 
dx = height/nHeight
length = nLength*dx
nLength_bnd = math.max(nHeight,nLength)+2
level = math.ceil(math.log(nLength_bnd)/math.log(2))
length_bnd = (2^level)*dx
dx_half = dx*0.5
dx_eps = length_bnd/2^20
zpos = dx_half
y_offset = -dx/2.0--height/2.0

-- Simulation parameters
--rho0_p = 1e3 --kg/m^3
---- density
--rho0_l = 1.0 
----kinematic viscosity
--nu_phy = 1e-6 --m^2/s
----dynamic viscosity
--mu = nu_phy*rho0_p --Pa s
-- Properties of sea water at 20C and 35g/kg salinity
rho0_p = 1025 --kg/m^3
mu = 1.08e-3 --Pa s
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
-- specific charge
charge_nr_H2O = 0.0
charge_nr_Na = 1.0
charge_nr_Cl = -1.0

ref_pot = -5e-3
pot_north = ref_pot
pot_south = ref_pot 
charge_north = 2e-5
charge_south = -2e-5
permit = 8.854e-12*80
conc_NaCl = 50 -- Molarity * 1e3
valence_sqr = 1
charge = 1.60217657e-19
k_b = 1.3805e-23
temp = 293
N_A = 6.02e23 
gasConst = 8.3144621
faraday = 96485.3365

-- Parameter for analytical solution
k = math.sqrt(2*conc_NaCl*valence_sqr*faraday^2/(permit*gasConst*temp))*height
k_H = k/height
print('k ', k, k_H)
--print(2e-5*math.cosh(k_H*height)/(permit*k_H*math.sinh(k_H*height)))
--print(2.0*permit*gasConst*temp/faraday*k_H*math.sinh(faraday*0/(2.0*gasConst*temp)))
function analy_pot_exp(x,y,z,t)  
  -- reciprocal of debye length
  ek = math.exp(k)
  eminusk = math.exp(-k)
  term_1 = (ek - 1)/(ek-eminusk)*math.exp(-k*y/height)  
  term_2 = (1-eminusk)/(ek-eminusk)*math.exp(k*y/height)
  if y>0 then
    return (term_1+term_2)*ref_pot 
  else
    return 0
  end  
end

local M = {}

local exp = math.exp

function M.cosh (x)
  if x == 0.0 then return 1.0 end
  if x < 0.0 then x = -x end
  x = exp(x)
  x = x / 2.0 + 0.5 / x
  return x
end


function analy_pot(x,y,z,t)
  k = math.sqrt(2*conc_NaCl*valence_sqr*faraday^2/(permit*gasConst*temp))*(height)
  return M.cosh(k*(y/height-0.5))/M.cosh(k/2)*ref_pot
end

function analy_pot_surfcharge(x,y,z,t)
  return charge_north/(permit*k_H)*M.cosh(k_H*y)/M.cosh(k/2.0)
end

function analy_electric_field(x,y,z,t)
  ek = math.exp(k)
  eminusk = math.exp(-k)
  term_1 = (ek - 1)/(ek-eminusk)*(-k)*math.exp(-k*y/height)  
  term_2 = (1-eminusk)/(ek-eminusk)*k*math.exp(k*y/height)
  return {0.0,-(term_1+term_2)*ref_pot/height,0.0}
end
--function analy_pot_2(x,y,z,t)
--  return charge_north*math.cosh(k_H*y)/(permit*k_H*math.sinh(k))
--end

--concentration of sea water
--salinity of sea water = 35
--http://en.wikipedia.org/wiki/Seawater
conc_Na = conc_NaCl --mol/m^3
conc_Cl = conc_NaCl --mol/m^3
rho_H2O = rho0_p - conc_NaCl*(mNa+mCl)
conc_H2O = rho_H2O/mH2O --mol/m^3
moleDens0 = conc_H2O+conc_Na+conc_Cl
---- mole fraction
moleFrac_H2O = conc_H2O/moleDens0
moleFrac_Na = conc_Na/moleDens0
moleFrac_Cl = conc_Cl/moleDens0
print('dx =', dx)

function analy_charge_density(x,y,z,t)
  pot = analy_pot(x,y,z,0)
  na = charge_nr_Na*conc_Na*math.exp(-faraday*charge_nr_Na/(gasConst*temp)*pot)
  cl = charge_nr_Cl*conc_Cl*math.exp(-faraday*charge_nr_Cl/(gasConst*temp)*pot)
  return (na+cl)*faraday
end
--print('charge_dens ', analy_charge_density(0,height/2.0,0,0), analy_pot(0,height/2.0,0))
