-- Musubi configuration file. 
package.path = package.path .. ";"..'$!seeder_path!$'.."/seeder.lua"
require 'seeder'

--flow parameters
Re = $!Re!$
if Re == 20 then
  u_in_phy = 0.3 --m/s
  steady_state = true
elseif Re == 100 then
  u_in_phy = 1.5 --m/s
  steady_state = false
else
  print('Reynolds number undefined')
end

u_mean_phy = 2.0*u_in_phy/3.0
nu_phy = 1e-3 --m^2/s
rho0_p = 1.0

nL = math.ceil(Dia/dx)
Re_check = u_mean_phy*Dia/nu_phy
-- set true for acoustic scaling, false for diffusive scaling
acoustic_scaling = false
rho0_l = 1.0
cs2_l = 1./3.
p0_l = rho0_l*cs2_l
if acoustic_scaling == true then
--acoustic scaling
  u_in_L = 0.05
  dt = u_in_L*dx/u_in_phy
  u_mean_L = 2.0*u_in_L/3.0
  nu_L = nu_phy*dt/dx^2
  omega = 1.0/(3.0*nu_L+0.5)
else
--diffusive scaling
  omega = 1.9
  nu_L = (1.0/omega-0.5)/3.0
  dt = nu_L*dx^2/nu_phy
  u_in_L = u_in_phy*dt/dx
  u_mean_L = 2.0*u_in_L/3.0
end

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

press_p = rho0_p*dx^2/dt^2
p0_p = press_p--*cs2_l

inlet = 'ubb'
outlet = 'expol'
tracking_fol = 'tracking/'
restart_fol = 'restart/'

-- This is a LUA script.
--dx = getdxFromLevel( {len_bnd=length_bnd, level=level})
--dt = getdtFromVel( {dx = dx, u_p = u_in_phy, u_l = u_in_L } )
--omega = getOmegaFromdt( {dx=dx, dt=dt, nu_p = nu_phy } )

-- Simulation name
simulation_name = 'channel2D'
mesh = '$!MESH!$'  -- Mesh information
printRuntimeInfo = false
scaling = 'acoustic'
control_routine = 'fast'

-- Time step settigs
tmax_p = 100 --s
tmax     =  tmax_p/dt --nSteps     -- total iteration number
interval = tmax/10
sim_control = {
  time_control = {
    max = { --tmax_p
       sim = tmax_p,
       clock = 2*60*60-5*60
    },
    interval = {iter=interval}
  },
  abort_criteria = {
    steady_state = steady_state
  }
}  

-- restart 
restart = {
      ead = restart_fol..'channel2D_lastHeader.lua',
      write = restart_fol,
      time_control = { min = tmax_p, max = {iter=tmax}, interval = {iter=tmax}}
 }


ebug = { debugMode = true, debugFiles = true,
          debugMesh = './debug/mesh_', debugStates = { 
  write = {
    folder    = './debug/',    -- the folder the restart files are written to
    interval  = 1,           -- dump restart file interval
    tmin      = 1,           -- first timestep to output
    tmax      = tmax+1       -- last timestep to output
    }
 }} 


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

fluid = { omega = omega, 
          rho0 = rho0_p }

interpolation_method = 'linear'                     

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

identify = {label='2D',layout='d3q19',kind='lbm_incomp'}
-- Boundary conditions
boundary_condition = {  
{ label = 'west', 
  kind = 'inlet_'..inlet, 
  --velocityX = u_in_phy,
  velocityX = { predefined = 'combined',
                transient= {predefined='smooth', min_factor = 0.25, 
                            max_factor=1.0, from_time=0, to_time=1000*dt}, 
                spatial = u_inflow
                                    
	      },
   velocityY = 0.0, velocityZ = 0.0},
{ label = 'east',
   kind = 'outlet_'..outlet,
   pressure = p0_p }, 
{ label = 'north', 
   kind = 'wall' },
{ label = 'south', 
   kind = 'wall' },
{ label = 'obst', 
   kind = 'wall_linearInterpolation', }
 }

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

-- lua variables
variable = {
  { name = 'coeff_fac', 
    ncomponents = 3, 
    vartype = 'st_fun',
    st_fun = {cD, cL, 0} 
  },
  { name  = 'coeff',
    ncomponents = 3, 
    vartype = 'operation',
    operation = { 
      kind = 'multiplication',
      input_varname = { 'bnd_force_phy', 'coeff_fac'}
    }  
  }
}

-- Tracking              
tracking = {
-- global flow
{ 
  label = 'global',
  folder = tracking_fol,
  variable = {
    'pressure_phy',
    'velocity_phy',
  },
  shape = { kind = 'global' },
  time_control = { min = tmax_p, max = tmax_p, interval = tmax_p },
  format = 'harvester'
},
-- pressure drop at two points over time: P(0.15,0.2,t) and P(0.25,0.2,t)
{
  label = 'probe_pressure',
  folder = tracking_fol,
  variable = {'pressure_phy'}, 
  shape = {
          {kind = 'canoND', object = {origin ={0.15-dx,0.2,zpos} }},
          {kind = 'canoND', object = {origin ={0.25+dx,0.2,zpos} }}
          },
  time_control = {min = 0, max = tmax_p, interval = 10*dt},
  format = 'ascii'      -- over time
},
-- convergence
 {
   label = 'convergence', 
   folder = tracking_fol,
   variable = {'pressure_phy'}, 
   shape = {
           {kind = 'canoND', object = {origin ={0.15-dx,0.2,zpos} }},
           {kind = 'canoND', object = {origin ={0.25+dx,0.2,zpos} }}
           },
   time_control = {min = 0, max = tmax_p, interval = 10*dt},
   format='convergence',
   reduction = 'average',
   convergence = {norm='average', nvals = 50, absolute = true,
   condition = { threshold = 1.e-10, operator = '<=' }}
 },
-- forces at the obstacle
-- {
-- label = 'force', 
-- folder = tracking_fol,
-- variable = {'bnd_force_phy','bnd_force'}, 
-- shape = { kind = 'boundary', boundary = {'obst'}},
-- time_control = {min = tmax, max = tmax, interval = tmax},
-- reduction = {'sum','sum'},
-- format = 'ascii'      --},
-- {
--  label = 'hvs_force', 
--  folder = tracking_fol,
--  variable = {'bnd_force_phy'}, 
--  shape = { kind = 'boundary', boundary = {'obst'}},
--  time_control = {min = tmax, max = tmax, interval = tmax},
--  format = 'harvester'      
-- },
-- velocity over time
{
 label = 'probe_velocity', 
 folder = tracking_fol,
 variable = {'velocity_phy'}, 
 shape = {
           {kind = 'canoND', object = {origin = {length*0.5,0.2,zpos}}}},
  time_control = {min = 0, max = tmax_p, interval = 10*dt},
  format = 'ascii'      
},
-- velocity over height
-- {
--   label = 'line_velocity_in', 
--   folder = tracking_fol,
--   variable = {'velocity_phy'}, 
--   shape = {kind = 'canoND', object = {origin = {dx_half,0.0,zpos},
--                                      vec = {0.0,height,0.0},
--                                      segments = nHeight+2}},
--   time_control = {min = tmax, max = tmax, interval = tmax},
--   format = 'asciiSpatial'      
-- },
-- {
--   label = 'line_velocity_inmean', 
--   folder = tracking_fol,
--   variable = {'velocity_phy'}, 
--   shape = {kind = 'canoND', object = {origin = {0.1/2.0,0.0,zpos},
--                                      vec = {0.0,height,0.0},
--                                      segments = nHeight+2}},
--   reduction = 'average',                                   
--   time_control = {min = 0, max = tmax, interval = tmax},
--   format = 'asciiSpatial'      
-- },
-- 
-- {
--   label = 'line_velocity_middle', 
--   folder = tracking_fol,
--   variable = {'velocity_phy'}, 
--   shape = {kind = 'canoND', object = {origin = {length*0.5,0.0,zpos},
--                                      vec = {0.0,height,0.0},
--                                      segments = nHeight+2 }},
--   time_control = {min = tmax, max = tmax, interval = tmax},
--   format = 'asciiSpatial'      
-- },
-- velocity over length
{
  label = 'line', 
  folder = tracking_fol,
  variable = {'pressure_phy','velocity_phy'}, 
  shape = {kind = 'canoND', object = {origin = {0.0,height/2.,zpos},
                                     vec = { length, 0.0,0.0},
                                     segments = nLength+2}
                                         },
  time_control = {min = tmax_p, max = tmax_p, interval = tmax_p},
  format = 'asciiSpatial'      
}

}
-- drag and lift coefficient at last time step
if Re == 20 then
  table.insert( tracking, {
    label = 'cyl_coeff',
    folder = tracking_fol,
    variable = {'coeff'},
    shape = { kind = 'boundary', boundary = {'obst'}},
    time_control = {min = tmax_p, max = tmax_p, interval = tmax_p},
    format = 'ascii',
    reduction = {'sum'}})
-- drag and lift coefficient over time    
elseif Re == 100 then
  table.insert( tracking,
  {
    label = 'cyl_coeff',
    folder = tracking_fol,
    variable = {'coeff'},
    shape = { kind = 'boundary', boundary = {'obst'}},
    time_control = { min = 0.0, max = tmax_p, interval = {iter=100} },
    format = 'ascii',
    reduction = {'sum'},
  --  transient_reduction = {'max'}
  })
end



