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

tMax = os.getenv("TMAX")
if tMax == nil then
  tMax = 256
end

nB = 3
balance = {
  dynamic = true,
  time_control = {
    min = { iter = tMax/(nB+1) },
    max = { iter = tMax - tMax/(nB+1) },
    interval = { iter = tMax/(nB+1) },
    check_iter = 64,
  },
}
dump_bc_timing = true

-- Simulation setup
logging = { level = 3 }
estart = {
  write = './restart/',
  time_control = {
    min = { iter = tMax },
    max = { iter = tMax },
    interval = { iter = tMax },
    -- check_iter = 64,
  },
}

--------------------------------------------------------------------------
Re = 10
scaling    = 'acoustic'
intp       = 'quadratic'
relaxation = 'bgk'
layout     = 'd3q19'
kind       = 'lbm'
inlet_kind  = 'inlet_ubb'
outlet_kind = 'outlet_expol_slow'
qValue = true
simulation_name = 'L'..level
folder = 'mesh_L'..level..'/'
mesh   = folder
flood_diagonal = false
printRuntimeInfo = false
smoothbounds = true
smoothlevels = true
use_sparse_alltoall = true

-- Geometry parameters
D = 2.0
R = D / 2
refine_level = 0
bounding_length = 64 * D
minLevel = level
maxLevel = level + refine_level
dxMax    = bounding_length / (2^minLevel)
dxMin    = bounding_length / (2^maxLevel)
dx       = dxMax
dx_half  = dxMin / 2
zpos     = dx_half
bounding_x = -15 * D - dxMax

L        = 40 * D
inlet_x  = bounding_x + dxMax/2
outlet_x = bounding_x + L - dxMax/2
top_y    = 15 * D + dxMax/2
bottom_y = -top_y
front_z  =  top_y
back_z   = -top_y
length   = outlet_x - inlet_x
H        = top_y - bottom_y
--------------------------------------------------------------------------

-- fluid property
rhoPhy = 1.0      -- kg / m^3
nuPhy  = 1.0    -- m^2 / s

-- simulaton parameters
-- Re = V * D / nu
U = Re * nuPhy / D
Ustart = 0.5*U
-- 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
A = math.pi * D * D / 4.0
cD = 2 / ( rhoPhy * U * U * A )
cL = 2 / ( rhoPhy * U * U * A )


cs2 = 1.0 / 3.0
cs  = math.sqrt(cs2)

omega = 1.95
nuLB  = ( 2.0/omega - 1.0 ) / 6.0
dt    = nuLB / nuPhy * dxMax * dxMax
ULB   = U * dt / dx
Mach  = ULB / cs

p0 = rhoPhy * dxMax^2 / dt^2 * cs2
p0 = 0
pOutlet = p0
