-- all units are in physical units
-- length (m)
-- mesh info
refinement_factor = 2
length = 2*math.pi
nLength = 60--refinement_factor*100
nWidth = 1
dx = length/nLength
nLength_bnd = nLength + 2
level = math.ceil(math.log(nLength_bnd)/math.log(2))
length_bnd = (2^level)*dx

-- simulation params
tmax_p = 0.05
tmax = 1000
-- time step size 
dt = tmax_p/tmax

-- molecular weights(kg/mol)
mH2O = 18.01528e-3 -- H20
mNa = 22.98977e-3 -- Na
mCl = 35.4527e-3  -- Cl
m_min = math.min(mH2O,mNa,mCl)

-- fake diffusivities
diff_diag = 1.0e-2
diff_H2O_Na = 1.31570699e-3
diff_H2O_Cl = 2.097388e-3
diff_Na_Cl = 2.95407e-5

-- density
rho0 = 1.0 --1000.0 -- kg/m^3
moleDens0 = 1.0

-- shear viscosity
nu_p = 1e-6 --m^2/s
--------------------------------------------------------
-- Initial condition for taylor green vortex testcase
--------------------------------------------------------
-- constant density for all species
-- kinematic pressure
function TGV_kinematicPressure(x,y,z)
  return (rho0/4.0) * (math.cos(4*math.pi*x/length) + math.cos(4*math.pi*y/length))
end

-- all species are initialized with mixture velocity
-- mixture velocityX
function TGV_velocityX(x,y,z)
  return math.sin(2*math.pi*x/length)*math.cos(2*math.pi*y/length)
end

-- mixture velocityY
function TGV_velocityY(x,y,z)
  return -math.cos(2*math.pi*x/length)*math.sin(2*math.pi*y/length)
end


-- molefraction initial conditions
function molefrac_H2O(x,y,z)
  --return 0.1*math.sin(2*math.pi*x/length)*math.cos(2*math.pi*y/length) + 0.3
  return 1.0/3.0
end

function molefrac_Na(x,y,z)
  --return 0.1*math.cos(2*math.pi*x/length)*math.sin(2*math.pi*y/length) + 0.3
  return 1.0/3.0
end

function molefrac_Cl(x,y,z)
  r1 = molefrac_H2O(x,y,z)
  r2 = molefrac_Na(x,y,z)
  return 1.0 - r1 - r2
end

-- velocity initial conditions for each species
-- Species 1 H2O
p1 = length/8 -- [0, L] 
function velX_H2O(x,y,z)
  return math.sin(2*math.pi*(x+p1)/length)*math.cos(2*math.pi*(y+p1)/length)
end
function velY_H2O(x,y,z)
  return -math.cos(2*math.pi*(x+p1)/length)*math.sin(2*math.pi*(y+p1)/length)
end

-- Species 2 Na
p2 = 0 -- 0 or L
function velX_Na(x,y,z)
  return math.sin(2*math.pi*(x+p2)/length)*math.cos(2*math.pi*(y+p2)/length)
end
function velY_Na(x,y,z)
  return -math.cos(2*math.pi*(x+p2)/length)*math.sin(2*math.pi*(y+p2)/length)
end

-- Species 3 Cl
function velX_Cl(x,y,z)
  chi1 = molefrac_H2O(x,y,z)
  chi2 = molefrac_Na(x,y,z)
  chi3 = 1.0 - chi1 - chi2
  r1 = chi1*moleDens0*mH2O
  r2 = chi2*moleDens0*mNa
  r3 = chi3*moleDens0*mCl
  rtotal = r1+r2+r3
  vX = TGV_velocityX(x,y,z)
  vX1 = velX_H2O(x,y,z)
  vX2 = velX_Na(x,y,z)
  return (rtotal/r3)*vX - ( (r1/r3)*vX1 + (r2/r3)*vX2 )
end

function velY_Cl(x,y,z)
  chi1 = molefrac_H2O(x,y,z)
  chi2 = molefrac_Na(x,y,z)
  chi3 = 1.0 - chi1 - chi2
  r1 = chi1*moleDens0*mH2O
  r2 = chi2*moleDens0*mNa
  r3 = chi3*moleDens0*mCl
  rtotal = r1+r2+r3
  vY = TGV_velocityY(x,y,z)
  vY1 = velY_H2O(x,y,z)
  vY2 = velY_Na(x,y,z)
  return (rtotal/r3)*vY - ( (r1/r3)*vY1 + (r2/r3)*vY2 )
end
