-- all units are in physical units
-- length (m)
-- mesh info
refinement_factor = 2
length = 2*math.pi
nLength = 60--refinement_factor*100
nWidth = 2
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.02
tmax = 800
-- 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)

-- diffusivities
-- 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
diff_diag = 1.0e-2
diff_H2O_Na = 1.31570699e-3
diff_H2O_Cl = 2.097388e-3
diff_Na_Cl = 2.95407e-3

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

-- shear viscosity
nu_p = 1e-6 --m^2/s
-- two smallness parameters to avoid division by zero
deltaRho1 = 0.001
deltaRho2 = 0.001
deltaRho3 = 0.001

---------------------------------------
-- Initial condition for dilute mixture
---------------------------------------
-- small parameter to control the behavior of the dominating solvent
tau_1 = 0.0001

function diluteH2O(x,y,z)
  return tau_1*(math.sin(2*math.pi*x/length) + deltaRho1 + 1.0 )
end

function diluteNa(x,y,z)
  p1 = tau_1*(math.sin(2*math.pi*x/length) + deltaRho1 + 1.0 )
  return 0.5 + deltaRho2 - 0.5 * p1
end

function diluteCl(x,y,z)
  p1 = tau_1*(math.sin(2*math.pi*x/length) + deltaRho1 + 1.0 )
  return 0.5 - deltaRho2 - 0.5 * p1
end

---------------------------------------
-- Initial condition for solvent mixture
---------------------------------------
-- small parameter to control the behavior of the dominating solvent
tau_2 = 0.001

-- phase shift of species 2 with respect to species 3
omega = 0.95*math.pi

function solventH2O(x,y,z)
  p1 = tau_2*(math.sin(2*math.pi*x/length + omega) + deltaRho2 + 1.0)
  p2 = tau_2*(math.sin(2*math.pi*x/length) + deltaRho3 + 1.0 )
  return 1.0 - p1 - p2
end

function solventNa(x,y,z)
  return tau_2*(math.sin(2*math.pi*x/length + omega) + deltaRho2 + 1.0)
end

function solventCl(x,y,z)
  return tau_2*(math.sin(2*math.pi*x/length) + deltaRho3 + 1.0 )
end


--------------------------------------------------------
-- Initial condition for maxwellian stefan tube testcase
--------------------------------------------------------
-- width of transition zone (width of transition from 0 to 1 is smeared across
-- two cells in x direction).
deltaX = length/20--10*length/nLength

-- constants for the first and second species' density profile
C1 = 0.1
C2 = C1

-- two smallness parameters to avoid divisons by zero 
deltaRho_st = (1.-C1-C2)/3.;

function stefanTubeH2O(x,y,z)
  if(x < length/2) then
    return C1*0.5*(1.0 - math.tanh((x-length/4)/deltaX)) + deltaRho_st
  else
    return C1*0.5*(1.0 + math.tanh((x-length/4-length/2)/deltaX)) + deltaRho_st
  end
end

function stefanTubeNa(x,y,z)
  if(x < length/2) then
    return C1*0.5*(1.0 + math.tanh((x-length/4)/deltaX)) + deltaRho_st
  else
    return C1*0.5*(1.0 - math.tanh((x-length/4-length/2)/deltaX)) + deltaRho_st
  end
end

function stefanTubeCl(x,y,z)
  if(x < length/2) then
    p1 = C1*0.5*(1.0 - math.tanh((x-length/4)/deltaX)) + deltaRho_st
    p2 = C1*0.5*(1.0 + math.tanh((x-length/4)/deltaX)) + deltaRho_st
    return 1.0 - p1 - p2
  else
    p1 = C1*0.5*(1.0 + math.tanh((x-length/4-length/2)/deltaX)) + deltaRho_st
    p2 = C1*0.5*(1.0 - math.tanh((x-length/4-length/2)/deltaX)) + deltaRho_st
    return 1.0 - p1 - p2
  end
end


--------------------------------------------------------
-- Initial condition for taylor green vortex testcase
--------------------------------------------------------
-- constant density for all species
TGV_density = 1.0/3.0

-- kinematic pressure
function TGV_kinematicPressure(x,y,z)
  return (rho0_p/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 condition
molefrac_H2O = solventH2O--TGV_density
molefrac_Na = solventNa--TGV_density
molefrac_Cl = solventCl--TGV_density

-- boundary condition
bc = {   
--                          { kind = 'wall',
--                            label = 'east'},
--                          { kind = 'wall',
--                            label = 'west'},
}
