require 'seeder'
physicalUnits = 'phyUnits'
--physicalUnits = 'LBUnits'
simulation_name = 'cylinder_level'..maxlevel..'_'..physicalUnits
mesh = 'mesh/'

tmax = 80000 + 1
time = {useIterations = true, max = tmax, interval = tmax/10, min = 1}

restart = {
  --read = 'restart/'..simulation_name..'_lastHeader.lua',
  --write = 'restart/',
  time = { min = tmax, max = tmax, interval = 1000},
  buffer = 8000000
}

--  PHysical units
densityPhy = 1000
viscPhy = 0.001
nuPhy = viscPhy / densityPhy
dpPhy = 0.01 -- pressure drop
lengthPhy = 0.08
diameterPhy = 0.02
radiusPhy = diameterPhy / 2
velPhy_mean = dpPhy * radiusPhy * radiusPhy / 8 / viscPhy / lengthPhy
velPhy_max = 2 * velPhy_mean
RePhy = velPhy_mean * diameterPhy / nuPhy

scaling = 'diffusive'
-- Lattice Units
omega = 1.8
cs2 = 1/3
rho0LB = 1
nuLB = 1.0/6.0*(2.0/omega-1.0)
dx = BBL/(2^maxlevel)
--dt_phy = getdtFromOmega( { dx = dx_phy, nu_p = nu_phy, omega = omega} )
dt = nuLB/nuPhy*dx*dx
velLB_mean = velPhy_mean * dt / dx
velLB_max = velLB_mean * 2
lengthLB = lengthPhy / dx
diameterLB = diameterPhy / dx
radiusLB = diameterLB / 2
dpLB = 8 * nuLB * velLB_mean * lengthLB / radiusLB / radiusLB
ReLB = velLB_mean * diameterLB / nuLB


print(' Use '..physicalUnits)
print('           dx = '..dx )
print('           dt = '..dt )
if physicalUnits == 'phyUnits' then
  physics = { dt = dt, rho0 = densityPhy, setPhysicalQuantities = true, }
  pressureInit = rho0LB * cs2 * densityPhy * dx * dx / dt / dt
  pressureInlet = pressureInit + dpPhy
  pressureOutlet = pressureInit
print('        nuPhy = '..nuPhy )
print('  diameterPhy = '..diameterPhy )
print('        dpPhy = '..dpPhy )
print('pressureInlet  = '..pressureInlet )
print('pressureOutlet = '..pressureOutlet )
print('  velPhy_mean = '..velPhy_mean )
print('   velPhy_max = '..velPhy_max )
print('   Renold Phy = '..RePhy )
else
  physics = { dt = dt, rho0 = densityPhy, setPhysicalQuantities = false, }
  pressureInit = rho0LB * cs2
  pressureInlet = pressureInit + dpLB
  pressureOutlet = pressureInit
print('    Renold LB = '..ReLB )
print('         nuLB = '..nuLB )
print('   diameterLB = '..diameterLB )
print('         dpLB = '..dpLB )
print('   velLB_mean = '..velLB_mean )
print('    velLB_max = '..velLB_max )
end

umax = 0.07914366
--Function definition for poiseuille profile across circular cross section
function circ_poise(x,y,z, t)
  return umax*(1 - ((x^2.0+y^2.0)/R^2.0))
end

identify = { label = physicalUnits, kind = 'lbm', relaxation = 'bgk', layout = 'd3q19' }
fluid = { omega = omega, rho0 = 1.0 }

initial_condition = {
  pressure = pressureInit,
  velocityX = 0.0,
  velocityY = 0.0,
  velocityZ = 0.0 } 

boundary_condition = {  
  { label = 'wall', kind = 'wall'},
--  { label = 'inlet', 
--    kind = 'inlet_ubb', 
--    velocityX = 0.0,
--    velocityY = 0.0,
--    velocityZ = {kind = 'combined', 
--                 transient = {predefined='smooth', min_factor=0.0, max_factor=1.0, from_time=0, to_time=tmax/10},
--                 spatial = circ_poise
--                },
--    density = 1.0
--  } ,  
  { label = 'inlet', 
    kind  = 'outlet_expol',
    pressure = pressureInlet},
  { label = 'outlet', 
    kind  = 'outlet_expol',
    pressure = pressureOutlet }
}

  -- Tracking
  tracking = {
    {
      label = 'global',
      variable = {
        'velocity_phy','pressure_phy','wss_phy','shearstress_phy',
        'velocity','pressure','wss','shearstress'
      },
      shape = {
       kind = 'global',
       object = {
         origin = { 0.01, 0.01, -4},
         vec = { { 0, 0, 8.0} },
         segments = { 300 }
       },
     },
      time = {min = tmax, max = tmax, interval = interval},
      format = 'harvester',
      folder = 'tracking/',
    },
    {
     label = 'velocity',
     variable = { 'velocity_phy' },
     shape = {
       kind = 'canoND',
       object = {
         center = {0.0, 0.0, -3.78},
         halfvec = { {0.99, 0.0, 0.0}, {0.0, 0.99, 0.0} },
         segments = {30, 30}
       },
     },
--     reduction = 'average',
     time = { min = tmax, max = tmax, interval = interval},
     format = 'asciiSpatial',
     folder = 'tracking/' 
    },
  } -- tracking table

