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

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

scaling = os.getenv("SCALING")
if scaling == nil then
  scaling = 'acoustic'
end

intp = os.getenv("INTP")
if intp == nil then
  intp = 'quadratic'
end

if os.getenv('SM') == "single_level" or os.getenv('SM') == nil then
  refine_level = 0
else
  refine_level = 1
end

tEnd = tonumber(os.getenv('TEND'))
if tEnd == nil then
  tEnd = 10
end
-- ramping time for inlet
-- tRamp = math.min(1, tEnd/5)
tRamp = tEnd/10
balance = {
  active = true,
  folder = 'balance/', -- directory where to store tmp data
  kind = 'simple', -- based on predefined scalefactor for each level
                  -- based on scaling (acoustic or diffusive)
  --kind  = 'levelwise', -- weights based on level
  dynamic = false,
  -- time = {
  --   useIterations = true, -- Timings given in iterations? default: true
  --   max = tmax,           -- Maximal iteration to reach
  --   interval = tmax+1,
  -- }    -- Interval for checking density
}

--------------------------------------------------------------------------
-- Geometry parameters
D = 2.0
R = D / 2
large_domain = false
if large_domain then
  inlet_x  = -15 * D
  outlet_x =  49 * D
  top_y    =  15 * D
else
  inlet_x  = -5 * D
  outlet_x = 17 * D
  top_y    =  5 * D
end
bottom_y = -top_y
front_z  = top_y
back_z   = -top_y
length   = outlet_x - inlet_x
L        = length
H        = top_y - bottom_y
bounding_length = L * 1.1

refine_level = 3
r1 = 0.1*D
r2 = 0.4*D
r3 = r2 + D

minLevel     = level
maxLevel     = level + refine_level
dxMax   = bounding_length / (2^minLevel)
dxMin   = bounding_length / (2^maxLevel)
dx_half = dxMin / 2
zpos    = dx_half
--------------------------------------------------------------------------

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

-- simulaton parameters
qValue = true
simulation_name = 'over_sphere'
folder = 'mesh/'
mesh   = folder
flood_diagonal = false
printRuntimeInfo = false
inlet_kind  = 'inlet_eq'
outlet_kind = 'outlet_pab'
smoothbounds = true
smoothlevels = true

-- 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 )

-- Simulation setup
logging = { level = 10 }

relaxation = 'mrt'
layout     = 'd3q27'
kind       = 'lbm'

interpolation_method = {
  method            = intp,
  test              = false,
  testEach          = false,
  testFluids        = false,
  noIntpFromFiner   = false,
  noIntpFromCoarser = false,
}

-- if scaling == 'acoustic' then
--   -- keep velocity constant, compute omega from dt
--   uLB = 0.1
--   dt  = uLB * dxMax
--   nuLB = nuPhy * dt/dxMax/dxMax
--   omega = 1.0 / ( 3.0 * nuLB + 0.5)
-- else
  -- keep omage constant, compute dt from omage
omega = 1.98
nuLB  = ( 2.0/omega - 1.0 ) / 6.0
dt    = nuLB / nuPhy * dxMax * dxMax
dx    = (bounding_length/(2^maxLevel))
ULB   = U * dt / dx
print(ULB)
nIters= math.ceil( 1 / dt )
-- end

cs2 = 1.0 / 3.0
p0 = rhoPhy * dxMax^2 / dt^2 * cs2
if kind == 'lbm_incomp' then
  p0 = 0
end
pOutlet = p0
