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

--scaling = 'diffusive'
scaling = 'acoustic'
interpolation_method = {
 method = 'quadratic',
}

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/Dia --s

u_max_phy = 3.0*u_mean_phy/2.0

u_max_L = 0.05
dt = u_max_L * dx / u_max_phy
nu_L = nu_phy * dt / dx^2.
omega = 1.0 / (3.0 * nu_L + 0.5 )

print('omega ', omega)
print('u_L ', u_max_L)
print('dt ', dt)

u_mean_L = 2.0*u_max_L/3.0
cs = math.sqrt(1.0/3.0)
machNr = u_max_L / cs
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,t) 
  --return (0.5/(rho0_p*nu_phy))*press_grad*y*(height-y)
  return {4.0*u_max_phy*y*(height-y)/height^2.0,0.0,0.0}
end

-- draf and lift coeff factors
cD = 2 / (rho0_p * u_mean_phy * u_mean_phy * Dia * dx)
cL = cD

-- Simulation name
simulation_name = 'channel'
mesh = 'mesh/' -- Mesh information
printRuntimeInfo = false
timing_file = 'mus_timing.res'

-- Time step settigs
if Re == 20 then
  tmax_p = 10 --s
else
  tmax_p = 10 --s
end
interval_p = 0.1 -- s
--tmax  = {iter=1000} 
--interval = {iter=100}
tmax = 2.0
interval = 0.1
sim_control = {
  time_control = {
    max = tmax,
    interval = interval
    --max = {iter=200},
    --interval = {iter=100}
  },
  abort_criteria = {
    stop_file = 'stop',
    steady_state = false
  }
}

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

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

fluid = { 
         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 = 'west', kind = 'inlet_ubb', velocity = 'inlet_vel'},
{ label = 'west', kind = 'inlet_ubb', velocity = u_inflow},
--{ label = 'west', kind = 'outlet_expol', pressure = ref_press+1e-3},
{ label = 'east', kind = 'outlet_expol', pressure = ref_press}, 
{ label = 'north', kind = bc_wall , fac = 1.0 },
{ label = 'south', kind = bc_wall}, 
{ label = 'obst', kind = bc_cyl},
 }

-- source term
sponge_start = 3.0/4.0*length
sponge_end = length
omega_min = 1.80
function sponge(x,y,z,t)
  if x>sponge_start then
    return omega+(omega_min-omega)*(x-sponge_start)/(sponge_end-sponge_start)
  elseif x<offset_cyl then 
    return omega_min+(omega-omega_min)*x/offset_cyl
  else
    return omega
  end
end 

glob_source = {
  sponge_omega = 'sponge'
}
-- lua variables
variable = {
  {
    name = 'sponge',
    ncomponents = 1,
    vartype = 'st_fun',
    st_fun = {
      fun = sponge,
      --fun = omega_min,
  ---- Shape for anonymous function has no effect    
      NOshape = {kind='canoND',
        object = {origin = {sponge_start,0.0,zpos},
                  vec = {  {length/4.0,0.0,0.0},
                           {0.0,height,0.0}
                  },
                  segments = {math.ceil(length/4.0/dx)+2, nHeight+2,0.0}
        }
      }  
    }  
  }
}

-- Tracking
tracking = {
{ label = 'vtk',
  folder = 'tracking/',
  variable = {'pressure_phy','velocity_phy'}, 
  shape = {kind = 'all'},
  --time_control = {min = 0, max = tmax_p, interval = {iter=50}},
  time_control = {min = 0, max = tmax, interval = interval},
  output= {  format = 'vtk'},
},
}

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