-- definition of the multispecies porous medium test case
level = 10
scale=1e-6 --m

interpolation_method = 'quadratic'
qVals = true
flood_diagonal = true
simname = 'porous'
-- geometry definition

-- give the bounding box origin and the size of the medium to simulate
porous_origin_original = { -0.806*scale, -0.794*scale, -0.793*scale }
porous_size   = { 1.59*scale, 1.59*scale, 1.59*scale }
-- 
-- geometry definition
height = math.max( porous_size[2], porous_size[3])
inlet_ext = 0.1*height
outlet_ext = 0.1*height
top_ext = 0.0*height
bottom_ext = 0.0*height
height_Y = porous_size[2]+top_ext+bottom_ext
height_Z = porous_size[3]
height_max = math.max(height_Y, height_Z)
length = height+inlet_ext+outlet_ext

cube_length = math.max( length, height )
length_bnd = 2*cube_length
-- calculate dx
dx = length_bnd/2^level
dx_half = dx*0.5
dx_offset = 2*dx

simname = 'porous'

-- 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 [m^2/s]
-- 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-9
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 = 2.0
-- time step size 
dt = (resi_diag/resi_latt)*dx*dx
tmax_p = 10.0 --sec
tmax = tmax_p/dt

-- Temperature in vapor phase
temp = 300 --K
--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
-- reference moleDens0
moleDens0 = 1000.0
--reynolds number
Re = 0.002
u_in_phy = Re*nu_phy/height_max

--lattice parameters
omega_diff = 2.0
nu_L = nu_phy*dt/dx^2
omega_kine = 3.0*nu_L
u_in_L = u_in_phy*dt/dx

----------------------------------------------
-- Initial conditions and Boundary conditions
----------------------------------------------
delta_mole = 0.05

function IC_velocity(x,y,z)
  --return 16.0 * u_in_phy * y * z * (height_Y-y)*(height_Z-z) / (height_Y^2.0*height_Z^2)
  return 0.0
end

function BC_velocity(x,y,z,t)
  return 16.0 * u_in_phy * y * z * (height_Y-y)*(height_Z-z) / (height_Y^2.0*height_Z^2)
  --return IC_velocity(x,y,z)
end

-- initial conditions
function IC_moleFrac_Na(x,y,z)
  return delta_mole--+delta_mole*y/height
end

function IC_moleFrac_Cl(x,y,z)
  return delta_mole--+delta_mole*(1-y/height)
end

function IC_moleFrac_H2O(x,y,z)
  return 1.0 - IC_moleFrac_Na(x,y,z) - IC_moleFrac_Cl(x,y,z)
end

-- boundary conditions
function BC_inlet_moleFlux_H2O(x,y,z,t)
  return {IC_moleFrac_H2O(x,y,z)*BC_velocity(x,y,z,t)*moleDens0,0.0,0.0}
end

function BC_inlet_moleFlux_Na(x,y,z,t)
  return {IC_moleFrac_Na(x,y,z)*BC_velocity(x,y,z,t)*moleDens0,0.0,0.0}
end

function BC_inlet_moleFlux_Cl(x,y,z,t)
  return {IC_moleFrac_Cl(x,y,z)*BC_velocity(x,y,z,t)*moleDens0,0.0,0.0}
end

function BC_top_moleFrac_Na(x,y,z,t)
  return 5*delta_mole
  --return IC_moleFrac_Na(x,y,z)
end

function BC_top_moleFrac_Cl(x,y,z,t)
  --return delta_mole
  return IC_moleFrac_Cl(x,y,z)
end

function BC_top_moleFrac_H2O(x,y,z,t)
  return 1.0-BC_top_moleFrac_Na(x,y,z)-BC_top_moleFrac_Cl(x,y,z)
  --return IC_moleFrac_H2O(x,y,z)
end

function BC_bottom_moleFrac_Na(x,y,z,t)
  --return delta_mole
  return IC_moleFrac_Na(x,y,z)
end

function BC_bottom_moleFrac_Cl(x,y,z,t)
  return 5*delta_mole
  --return IC_moleFrac_Cl(x,y,z)
end

function BC_bottom_moleFrac_H2O(x,y,z,t)
  --return delta_mole
  --return IC_moleFrac_H2O(x,y,z)
  return 1.0-BC_bottom_moleFrac_Na(x,y,z)-BC_bottom_moleFrac_Cl(x,y,z)
end

-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
--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_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,t)*moleDens0 
end


function moleflux_Na(x,y,z,t)
  return IC_Na(x,y,z)*BC_velocity(x,y,z,t)*moleDens0
end
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
