-- Musubi configuration file. 
-- This is a LUA script.
require 'seeder'
identify = {label='3D',layout='d3q19', relaxation='bgk', kind = 'lbm_incomp' }

scaling = 'acoustic'
interpolation_method = {  method = 'weighted_average', }

logging = {level=5}

if qVal then
  bc_cyl = 'wall_linearInterpolation'
else
  bc_cyl = 'wall'
end
bc_wall = 'wall'

--flow parameters
Re = 100
rho0_p = 1.0
rho0_L = 1.0
cs2 = 1.0/3.0
nu_phy = 1e-3 --m^2/s
u_mean_phy = Re*nu_phy/height --s

u_max_phy = 3.0*u_mean_phy/2.0

if scaling == 'acoustic' then
  u_max_L = 0.0375
  dt = u_max_L * dx / u_max_phy
  nu_L = nu_phy * dt / dx^2.
  omega = 1.0 / (3.0 * nu_L + 0.5 )
elseif scaling == 'diffusive' then
  print 'ERROR: scaling not supported for multilevel'
  os.exit()
else
  print 'ERROR: scaling not supported'
  os.exit()
end

u_mean_L = 2.0*u_max_L/3.0
cs = math.sqrt(1.0/3.0)
machNr = u_max_L / cs
print( 'machNr =', machNr )
print( 'Omega =', omega )
if machNr > 0.1 then
  print ('ERROR: Mach number is above LBM limit. Ma =',machNr)
  os.exit()
end

ref_press = 1.0
press_drop = 8*u_max_phy*rho0_p*nu_phy*length/height^2
press_grad = press_drop/length

function u_inflow(x,y,z)
  return 4.0*u_max_phy*y*(height-y)/height^2.0
end

-- Simulation name
simulation_name = 'channel_D3Q19_avg_incomp'
mesh = 'mesh/' -- Mesh information
printRuntimeInfo = false

-- Time step settigs
if Re == 20 then
  tmax_p = 20 --s
else
  tmax_p = 50 --s
end
tmax  =  tmax_p/dt    -- total iteration number
interval_p = 1 -- s
sim_control = {
  time_control = {
    max = {sim=tmax_p},
    interval = interval_p,
  },
  abort_criteria = {
    stop_file = 'stop',
    steady_state = true,
    convergence = {
      variable = {'pressure_phy', 'velocity_phy'},
      shape = 'all',
      time_control = {min = 0, interval = {iter=10}},
      reduction = {'average', 'average'},
      norm='average', nvals = 10, absolute = true,
      condition = {
        {threshold = 1.e-5, operator = '<=' },
        {threshold = 1.e-5, operator = '<=' }
      }
    }
  }
}

nodebug= {debugMode = true, debugFiles=true, debugMesh='debug/',
          logging = {level=3, filename = 'dbg', root_only = true },
         }
-- logging

-- needed to dump variable in physical unit
physics = { dt = dt, rho0 = rho0_p }

fluid = {kinematic_viscosity=nu_phy,
         omega = omega, rho0 = rho0_p }

-- Initial condition
initial_condition = { pressure = ref_press,
                      velocityX = 0.0,
                      velocityY = 0.0,
                      velocityZ = 0.0, }


-- Boundary conditions
boundary_condition = {
{ label = 'north', kind = bc_wall , fac = 1.0 },
{ label = 'south', kind = bc_wall},
{ label = 'obst', kind = bc_cyl},
{ label = 'west', kind = 'inlet_ubb', velocity = 'inlet_vel'},
{ label = 'east', kind = 'outlet_expol', pressure = ref_press},
 }

-- lua variables
variable = {
  { name = 'inlet_vel',
    ncomponents = 3,
    vartype = 'st_fun',
    st_fun = {
      predefined = 'combined',
      temporal = {predefined='smooth', min_factor = 0,
                  max_factor=1.0, from_time=0, to_time=1},
      spatial = {
        predefined='parabol',
        shape = {
          kind = 'canoND',
          object = {
            origin  = {0.0,0.0,dx/2.0},
            vec = {0.0,height,0.0}
          },
        }, -- shape
        amplitude = {u_mean_phy,0.0,0.0}
      }
    }

  },
}

-- Tracking
tracking = {
--  { label = 'vtk',
--    folder = 'tracking/',
--    variable = {'pdf','density','density_phy','pressure_phy','pressure','velocity_phy','velocity','shear_stress','process'},
--    shape = {kind = 'all'},
--    time_control = { min = tmax_p, max = tmax_p },
--    output= {  format = 'vtk'},
--  },
  { label = 'ref_line',
    variable = { 'velocity', 'pressure' },
    folder = 'tracking/',
    shape = { kind = 'canoND', object = {origin = {0.05, 0.2, 0.005},
                                       vec = {1.0, 0.0, 0.0},
                                       segments = 11,
                                       distribution = 'equal'}},
    output = { format = 'asciiSpatial' },
    time_control = { min = tmax_p, interval = tmax_p }
  },
}
-- end of tracking table


norestart = {
      ead = 'restart/'..'/channel2D_lastHeader.lua',
      write = 'restart/',
      time_control = { min = tmax_p, max = tmax_p, interval = tmax_p}
 }
