--This file contains common parameters

level = tonumber(os.getenv('LEVEL'))
if level == nil then
  level = 6
end
Re = tonumber(os.getenv('Re'))
if Re == nil then
  Re = 30
end

Mach = tonumber(os.getenv('Mach'))
if Mach == nil then
  Mach = 0.1
end

if os.getenv('MULTI') == 'y' then
  multi_level = 1
else
  multi_level = 0
end

if os.getenv('SCALING') == 'acoustic' then
  diffusive = false
  scaling = 'acoustic'
else
  diffusive = true  -- diffusive or acoustic scaling?
  scaling = 'diffusive'
end

if os.getenv('INTP') == 'linear' then
  interpolation_method = 'linear'
else
  interpolation_method = 'quadratic'
end

identify_kind = os.getenv('KIND')
if identify_kind == nil then
  identify_kind = 'lbm'
end

simulation_name = 'channel2D_L'..level
relaxation = 'mrt'
useVelocityInlet = true
usePeriodic = true -- true = 2D channel, false = 3D channel
printRuntimeInfo = false
layout = 'd2q9'

mesh = 'mesh/' -- Mesh files location

length = 4.0
dx = length / 2.^level
height = length / 4
if useVelocityInlet then
  inlet_kind = 'inlet_ubb'
else
  inlet_kind = 'outlet_expol'
end
outlet_kind = 'outlet_expol'
wall_kind   = 'wall'
top_y = height/2 + dx/2
bottom_y = -height/2 - dx/2
ly = height / dx       -- number of elements for channel height
lx = length / dx - 2   -- number of elements for channel height
R  = height / 2
inlet_x  = -length/2 + dx/2
outlet_x =  length/2 - dx/2

-- Physical parameters
densityPhy = 1
viscPhy = 1
nuPhy = viscPhy / densityPhy

-- reference LB quantities
rho0 = 1.
cs2  = 1./3.
cs   = math.sqrt(cs2)
velMaxLB  = Mach * cs
nuLB = velMaxLB * ly / Re
omega = 1.0 / ( 3.0 * nuLB + 0.5 )
dt   = nuLB * dx * dx / nuPhy
velMaxPhy = velMaxLB * dx / dt
p0LB = rho0*cs2
tau  = 1.0 / omega
tmaxLB  =  2000 * ( 4.^(level - 5 ))
tmax = tmaxLB

p0Phy = p0LB * densityPhy * dx * dx / dt / dt -- pressure at outlet
p0Phy = 0
pOutlet = p0Phy
rhoLBOutlet = pOutlet * dt * dt / densityPhy / dx / dx / cs2
tmaxPhy = tmaxLB * dt

dp_dx = -2.0 * velMaxPhy * viscPhy / R / R
dp = dp_dx * ( length - 3*dx )
pInlet = pOutlet - dp

-- Analytical solutions
function wssA(x, y, z, t)
  return viscPhy * 2.0 * math.abs(y) * velMaxPhy / R / R
end
function strainRateA(x, y, z,t)
  return -1*(y) * velMaxPhy / R / R
end
function velA(x, y, z,t)
  return velMaxPhy * ( 1 - (y)^2 / R^2 )
end
-- outlet pressure is @ length/2 - dx
-- outlet STL is at length/2 - dx/2,
-- the outlet element is at length/2 - dx/2 - dx
function pressureOutlet(x,y,z)
  return pOutlet
end
-- pressure analytical
function pressureA(x,y,z,t)
  return pOutlet + dp_dx * ( x - (length/2-3*dx/2) )
end
