-- Musubi configuration file. 
require 'seeder'

--flow parameters
nu_phy = 5e-2 --m^2/s
rho0_p = 1.0
rho0_L = 1.0
cs2 = 1.0/3.0

Re = 10
u_mean_phy = Re*nu_phy/height
u_max_phy = 3.0*u_mean_phy/2.0

ref_press = 1.0
--press_drop = 1e-4
press_drop = 8*u_max_phy*rho0_p*nu_phy*(length)/height^2
press_grad = press_drop/(length)

function u_analy(x,y,z) 
  return (0.5/(rho0_p*nu_phy))*press_grad*y*(height-y)
  --return (0.5/(rho0_p*nu_phy))*press_grad*(height*0.5-y)*(height*0.5+y)
end

function u_inflow(x,y,z,t) 
  return {u_analy(x,y,z), 0.0, 0.0}
  --return {u_mean_phy, 0.0,0.0}
end

u_max_phy = u_analy(0.0,height/2.0,0.0)
u_mean_phy = 2.0*u_max_phy/3.0

Re = u_mean_phy*height/nu_phy

scaling = 'acoustic'

--scaling = 'diffusive'
--if scaling == 'acoustic' then
----acoustic scaling
--  u_max_L = 0.001
--  dt = u_max_L*dx/u_max_phy
--  u_mean_L = 2.0*u_max_L/3.0
--  nu_L = nu_phy*dt/dx^2.
--  omega = 1.0/(3.0*nu_L+0.5)
--else
--diffusive scaling
  omega = 0.9
  nu_L = (1.0/omega-0.5)/3.0
  dt = nu_L*dx^2/nu_phy
  u_mean_L = u_mean_phy*dt/dx
  u_max_L = 3.0*u_mean_L/2.0
--end

-- Analytical solutions
function wssA(x, y, z)
  return nu_phy * 2 * math.abs(y-height*0.5) * u_max_phy / (height*0.5)^2.0
end
function strainRateA(x, y, z)
  return -1*(y) * u_max_phy / R / R
end
function pressA(x,y,z,t)
  return ref_press+press_drop*((length)-x)/(length)
end

inletBC = 'inlet_ubb'--'press_dirichlet'
outletBC = 'outlet_expol'--'press_dirichlet'
curved = false
qVal_press = 0.0
qVal_vel = 0.5
wallBC = 'vel_dirichlet'
wallBC = 'wall'

-- Simulation name
simulation_name = 'channel'
mesh = 'mesh/' -- Mesh information
printRuntimeInfo = false
control_routine = 'fast'
NOdebug = {logging = {level=1, filename='dbg', root_only=false}}
-- Time step settigs
tmax_p = 10000 --s
tmax =  tmax_p/dt    -- total iteration number
interval_p = 1 -- s
interval = interval_p/dt
sim_control = {
  time_control = { 
    max = {sim=tmax_p}, 
    interval = {sim=interval_p}
   },
  abort_criteria = {
    stop_file = 'stop',
    steady_state = true,
    convergence = {
     variable = {'pressure_phy','vel_mag_phy'}, 
     shape = {kind = 'all'},
     time_control = {min = 1, max = tmax_p, interval = 10*dt},
     reduction = {'average','average'},
     norm='average', nvals = 50, absolute = true,
     condition = {
       { threshold = 1.e-10, operator = '<=' },
       { threshold = 1.e-10, operator = '<=' }}
    },
  }
}

NOrestart = {
      ead = 'restart/channel2D_lastHeader.lua',
      write = 'restart/',
}

fluid = { 
           omega = omega,
           rho0 = rho0_p }

interpolation_method = 'linear'                     

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

identify = {
  label='2D',
  layout='d2q9', 
  relaxation='mrt',
  kind = 'lbm_incomp'
}

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

-- Boundary conditions
boundary_condition = {  
  {
    label = 'west', 
    kind = inletBC,
    velocity = u_inflow,
    --pressure = ref_press+press_drop,
    pressure = pressA,
    curved = false,
    qVal = qVal_vel,
  },  
  {
    label = 'east',
    kind = outletBC,
    pressure = pressA,
    --pressure = ref_press,
    velocity = u_inflow,
    curved = false,
    qVal = qVal_press
  }, 
  {   
    label = 'north', 
    kind = wallBC,
    velocity = {0.0,0.0,0.0},
    qVal = qVal_vel, 
    curved = true,
  },
  {
    label = 'south', 
    kind = wallBC,
    velocity = {0.0,0.0,0.0},
    qVal = qVal_vel, 
    curved = true,
  },
}

variable = {
  { name='vel_an', 
    ncomponents=1, 
    vartype = 'st_fun', 
    st_fun=u_analy 
  },
  { name='wss_an', 
    ncomponents=1, 
    vartype = 'st_fun', 
    st_fun=wssA 
  },
  { name='press_an', 
    ncomponents=1, 
    vartype = 'st_fun', 
    st_fun=pressA 
  },
  { name='diff_press', 
    ncomponents=1, 
    vartype = 'operation', 
    operation = {
      kind='difference', 
      input_varname={'pressure_phy','press_an'} 
    }
  },
  { name='diff_vel', 
    ncomponents=1, 
    vartype = 'operation', 
    operation = {
      kind='difference', 
      input_varname={'vel_mag_phy','vel_an'} 
    }
  },
  { name='diff_wss', 
    ncomponents=1, 
    vartype = 'operation', 
    operation = {
      kind='difference', 
      input_varname={'wss_phy','wss_an'} 
    }
  },
}

-- Tracking              
tracking = {
--  {
--    label = 'probe', 
--    folder = 'tracking/',
--    variable = {'pressure_phy','velocity_phy'}, 
--    shape = {
--      kind = 'canoND', 
--      object = {
--        origin ={length*0.5,height*0.5,zpos} 
--      } 
--    },
--    time_control = {min= {iter= 1}, max = {iter= tmax}, interval = {iter= 10}},
--    output = {format = 'ascii'}
--  },
  {
    label = 'vtk', 
    folder = 'vtkfiles/',
    variable = {'pressure_phy','velocity_phy','press_an','vel_an'}, 
    shape = {kind = 'all'},
    time_control = {min= 0, max = tmax_p, interval = tmax_p/100},
    output = {format = 'vtk'}
  },
--  {
--    label = 'pressInlet', 
--    folder = 'tracking/',
--    variable = {'pressure_phy'}, 
--    shape = {
--      kind = 'canoND', 
--      object = {
--        origin ={dx/2.0,0.0,zpos},
--        vec = {0.0,height,0.0},
--        segments = nHeight+2
--      } 
--    },
--    time_control = {min = {iter= tmax}, max = {iter= tmax}, interval = {iter= tmax+1}},
--    output = {format = 'asciiSpatial'}
--  },
--  {
--    label = 'pressOutlet', 
--    folder = 'tracking/',
--    variable = {'pressure_phy'}, 
--    shape = {
--      kind = 'canoND', 
--      object = {
--        origin ={length-dx/2.0,0.0,zpos},
--        vec = {0.0,height,0.0},
--        segments = nHeight+2
--      } 
--    },
--    time_control = {min = {iter= tmax}, max = {iter= tmax}, interval = {iter= tmax+1}},
--    output = {format = 'asciiSpatial'}
--  },
--  {
--    label = 'press_line', 
--    folder = 'tracking/',
--    variable = {'pressure_phy','press_an','diff_press'}, 
--    shape = {
--      kind = 'canoND', 
--      object = {
--        origin ={0.0,height*0.5,zpos},
--        vec = {length,0.0,0.0},
--        segments = nLength+2
--      } 
--    },
--    time_control = {min = {iter= tmax}, max = {iter= tmax}, interval = {iter= tmax+1}},
--    output = {format = 'asciiSpatial'}
--  },
--  { label = 'press_l2norm',
--    variable = {'diff_press', 'press_an'},
--    shape = {
--      kind = 'canoND', 
--      object = {
--        origin ={0.0,height*0.5,zpos},
--        vec = {length,0.0,0.0},
--        segments = nLength+2
--      }
--    },
--    time_control = { min = {iter=tmax}, max = {iter=tmax}, interval = {iter=tmax+1} },
--    reduction = {'l2norm', 'l2norm'},
--    folder = './tracking/',
--    output = {format = 'asciiSpatial'}
--  },
--  {
--    label = 'vel_centerLine', 
--    folder = 'tracking/',
--    variable = {'velocity_phy'}, 
--    shape = {
--      kind = 'canoND', 
--      object = {
--        origin ={0.0,height*0.5,zpos},
--        vec = {length,0.0,0.0},
--        segments = nLength+2
--      } 
--    },
--    time_control = {min = {iter= tmax}, max = {iter= tmax}, interval = {iter= tmax+1}},
--    output = {format = 'asciiSpatial'}
--  },
--  { label = 'vel_spatial',
--    variable = {'velocity_phy','vel_mag_phy','vel_an', 'diff_vel'},
--    shape = {
--      kind = 'canoND', 
--      object = {
--        origin ={length*0.5,0.0,zpos},
--        vec = { 0.0, height, 0.0},
--        segments = nHeight+2
--      }
--    },
--    time_control = { min = {iter=tmax}, max = {iter=tmax}, interval = {iter=tmax+1} },
--    folder = './tracking/',
--    output = {format = 'asciiSpatial'}
--  },
--  { label = 'vel_l2norm',
--    variable = {'diff_vel','vel_an'},
--    shape = {
--      kind = 'canoND', 
--      object = {
--        origin ={length*0.5,0.0,zpos},
--        vec = { 0.0, height, 0.0},
--        segments = {nHeight+2} 
--      }
--    },
--    time_control = { min = {iter=tmax}, max = {iter=tmax}, interval = {iter=tmax+1} },
--    reduction = {'l2norm', 'l2norm'},
--    folder = './tracking/',
--    output = {format = 'asciiSpatial'}
--  },
--  { label = 'wss_spatial',
--    variable = { 'wss_phy', 'wss_an', 'diff_wss'},
--    shape = {
--      kind = 'canoND', 
--      object = {
--        origin ={length*0.5,0.0,zpos},
--        vec = { 0.0, height, 0.0},
--        segments = {nHeight+2} 
--      }
--    },
--    time_control = { min = {iter=tmax}, max = {iter=tmax}, interval = {iter=tmax+1} },
--    folder = './tracking/',
--    output = {format = 'asciiSpatial'}
--    },
--  
--  { label = 'wss_l2norm',
--    variable = {'diff_wss','wss_an'},
--    shape = {
--      kind = 'canoND', 
--      object = {
--        origin ={length*0.5,0.0,zpos},
--        vec = { 0.0, height, 0.0},
--        segments = {nHeight+2} 
--      }
--    },
--    time_control = { min = {iter=tmax}, max = {iter=tmax}, interval = {iter=tmax+1} },
--    reduction = {'l2norm', 'l2norm'},
--    folder = './tracking/',
--    output = {format = 'asciiSpatial'}
--  },
}
