level = os.getenv("LEVEL")
if level == nil then
  level = 10
end

Re = tonumber( os.getenv("RE") )
if Re == nil then
  Re = 20
end

tEnd = tonumber(os.getenv('TEND'))
if tEnd == nil then
  tEnd = 12
end

-- Geometry parameters
H   = 0.41
l_h = 5
L   = H * l_h
bounding_length = 2.56

D = 0.1
R = D / 2

dx       = bounding_length / (2^level)
dx_half  = dx / 2
zpos     = dx_half

cylinder_x = 0
cylinder_y = 0.0
top_y      =     4.2 * R + dx_half
bottom_y   = - ( 4.0 * R + dx_half )
inlet_x    = - ( 4.0 * R + dx_half )
outlet_x   = inlet_x + L
periodic_z1 = -dx_half
periodic_z2 =  dx_half + dx

-- fluid property
rhoPhy  = 1.0             -- kg / m^3
nuPhy   = 0.001           -- m^2 / s
viscPhy = nuPhy * rhoPhy  -- kg / ( m s )

-- simulaton parameters
qValue = true
if qValue then
  bc_wall = 'wall_libb'
else
  bc_wall = 'wall'
end
relaxation = 'mrt'
simulation_name = 'L'..level
folder = 'mesh_L'..level..'/'
-- folder = 'mesh/'
mesh   = folder
flood_diagonal = false
printRuntimeInfo = false
inlet_kind  = 'vel_dirichlet'
outlet_kind = 'press_dirichlet'
smoothbounds = true
smoothlevels = true
scaling = 'acoustic'

-- Re=20,  Um = 0.3
-- Re=100, Um = 1.5
Um = 0.3 * Re / 20
Uaverage = Um * 2 / 3
-- Lift or Drag coefficient are defined as
-- Cl = L / ( 0.5 * rho * U * U * A )
-- Cd = D / ( 0.5 * rho * U * U * A )
-- where L or D are force, U is velocity, A is area
-- from tracking we get L or D, then multiply with cD or cL, we can get Cl, Cd
cD = 2 / ( rhoPhy * Uaverage * Uaverage * D * dx )
cL = 2 / ( rhoPhy * Uaverage * Uaverage * D * dx )

function inlet_U( x, y, z, t )
  return 4.0 * Um * ( y-ref_y ) * ( H - (y-ref_y) ) / H / H
end

-- Simulation setup
logging = { level = 5 }
layout  = 'd2q9'
kind    = 'lbm'

cs2   = 1./3.
cs    = math.sqrt(cs2)
omega = 1.8
nuLB  = ( 1.0 / omega - 0.5 ) / 3.0
dt    = nuLB / nuPhy * dx * dx
uLB   = Um / dx * dt
Mach  = uLB / cs
velMaxLB  = Mach * cs
velMaxPhy = velMaxLB * dx / dt
tRamp = 1
p0 = 0
pOutlet = p0
