-- Musubi configuration file. 
-- ------------------------------
-- Steady flow in a round tube --
-- ------------------------------
--
--package.path = package.path .. ";"..'$!seeder_path!$'.."/seeder.lua"
require 'seeder'

-- parameters taken from common.lua

---------------------
-- flow parameters --
---------------------

nu_phy = 1e-3 --m^2/s
rho0_p = 1.0
rho0_L = 1.0
cs = math.sqrt(1.0/3.0)

ref_press = 1.0 				
press_drop = 0.1				
press_grad = press_drop/length
press_inlet = ref_press+press_drop

u_in_phy = radius^(2.0) / (4.0*nu_phy*rho0_p) * press_grad
u_mean_phy = 0.5*u_in_phy -- u_max

Re = u_mean_phy*2.0*radius / nu_phy
print('Re = ', Re)
print('u_in_phy = ', u_in_phy)


---------------------------
acoustic_scaling = false --	-- set true for acoustic scaling, false for diffusive scaling
---------------------------

if acoustic_scaling == true then
----------------------
-- acoustic scaling --
----------------------
  u_in_L = 0.08
  dt = u_in_L*dx/u_in_phy
  u_mean_L = 0.5*u_in_L
  nu_L = nu_phy*dt/dx^2.
  omega = 1.0/(3.0*nu_L+0.5)
else
-----------------------
-- diffusive scaling --
-----------------------

  omega = 1.95 --$!omega!$
  nu_L = (1.0/omega-0.5)/3.0
  --nu_L = 0.0500
  --omega = 1.0/(3.0*nu_L+0.5)
  dt = nu_L*dx^2/nu_phy
  u_mean_L = u_mean_phy*dt/dx
  u_in_L = 2.0*u_mean_L
end

cutoff_min = length/5.0
cutoff_max = length
-- max omega cut off ratio
w_max = 1.0
-- min omega cut off ratio
w_min = 10.0
function omega_spatial(x,y,z)
  if ( x < cutoff_min ) then
    res = w_max
  elseif (x >= cutoff_max)  then
    res = w_min
  else
    slope = (w_max-w_min)/(cutoff_min-cutoff_max)
    res =  x*slope + w_max - cutoff_min*slope
  end
  return res
end

--------------------------
-- Analytical solutions --
--------------------------

function u_analy(x,y,z,t)
  return radius^(2.0) / (4.0*nu_phy*rho0_p) * press_grad * (  1.0 - (y^(2.0)+z^(2.0)) / (radius^(2.0))  )
end

wss_max = 4.0*nu_phy*rho0_p*u_mean_phy / radius

function wssA(x, y, z,t) -- tau = mu*(du/dr)
  return nu_phy*rho0_p*( radius^(2.0) / (4.0*nu_phy*rho0_p) * press_grad * (2.0*math.sqrt(y^(2.0)+z^(2.0)) / (radius^(2.0))  ) )
end

--function strainRateA(x, y, z, t)
--  return -1*(y) * u_in_phy / R / R
--end

function pressA(x,y,z, t)
  return ref_press+press_drop*(length-x)/length
end

----------------
-- X velocity --
----------------

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

--            { predefined = 'combined',
--            transient= {predefined='smooth', min_factor = 0.5, max_factor=1.0, 
--            from_time=0, to_time=1}, 
--            spatial = u_mean_phy
--            }


folderRead = 'restart/'
folderWrite = 'restart/'

---------------------
-- Simulation name --
---------------------

simulation_name = 'pipe'
logging = {level=4}
mesh = 'mesh/'--'$!MESH!$' ---- Mesh information
timing_file = 'mus_timing.res'
printRuntimeInfo = false
control_routine = 'fast'
require_moments = true

-----------------------
-- Time step settigs --
-----------------------

tmax_p = 10 --s
tmax   = tmax_p/dt    -- total iteration number

interval_p = tmax_p/50
interval   = interval_p/dt
wall_clock = 4*60*60-10*60
sim_control = {
  time_control = { 
    max = {sim=tmax_p,iter=tmax,clock=wall_clock}, 
    interval = {iter=interval}
   },
   abort_criteria = {
     stop_file = 'stop',
     steady_state = true,
     convergence = {
       variable = {'pressure_phy'}, 
       shape = {kind = 'all'}, 
       time_control = {min = 0, max = tmax_p, interval = 10*dt},
       reduction = {'average'},
       norm='average', nvals = 100, absolute = true,
       condition = { threshold = 1.e-15, operator = '<=' }
     }
   }  
} --sim_control

----------------------
-- restart          --
----------------------

restart = {
  NOread = folderRead,
  write = folderWrite,
  time_control={min=tmax_p}
}

------------------------------------------------
--balance = { active = false, dynamic = true, kind = 'levelwise',
--            time_control = { interval = {iter=3}, min = {iter=14}, max = {iter=tmax}}
--  }

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

-----------------------
-- fluid information -- 
-----------------------

fluid = { 
        --omega_ramping = {predefined='smooth', min_factor = 0.5, 
        --                   max_factor=1.0, from_time=0, to_time=10000.0*dt}, 
           --kine_shear_viscosity = nu_phy,
           --visc_spatial = { fun = omega_spatial, store = true},
           bulk_viscosity = 0.01,
           omega = omega,
           rho0 = rho0_p }

interpolation_method = 'linear'                     

------------------------
-- Initial conditions --
------------------------

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

identify = {label='bgk',layout='d3q27', relaxation='mrt',
            kind = 'lbm'
}

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

----------------------
-- Inlet and Outlet --
----------------------

inletBC = 'inlet_ubb'
outletBC = 'outlet_expol'

-------------------------
-- Boundary conditions --
-------------------------

boundary_condition = {
  { 
    label = 'inlet',
    kind = inletBC,
    --pressure = press_inlet,
    velocity = u_inflow,
  },
  {
    label = 'outlet',
    kind = outletBC,
    pressure = ref_press
  }, 
  {
    label = 'cylinder',
    kind = 'wall_libb',
  }
}

---------------
-- variables --
---------------

variable = {
  --analytic solution
  { 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 },
    --error
  { name = 'diff_pressure', 
    ncomponents = 1, 
    vartype = 'operation',
    operation = {
      kind = 'difference',
      input_varname = {'pressure_phy','press_an'}
    }
  },
  { name = 'diff_wss', 
    ncomponents = 1, 
    vartype = 'operation',
    operation = {
      kind = 'difference',
      input_varname = {'wss_phy','wss_an'},
    }
  },
  { name = 'diff_vel', 
    ncomponents = 1,
    vartype = 'operation',
    operation = {
      kind = 'difference',
      input_varname = {'vel_mag_phy','vel_an'}
    }
  },
}--variable


--------------
-- Tracking --
--------------
              
tracking = {
{
  label = 'vtk', 
  folder = 'tracking/',
  variable = {'pressure_phy','velocity_phy','vel_an'}, 
  shape = {kind = 'all'},
  time_control = {min= {iter=0}, max = {iter= tmax}, interval = tmax_p/10},
  output = {format = 'vtk'}      
},
{
  label = 'probe_inlet', 
  folder = 'tracking/',
  variable = {'pressure_phy','velocity_phy'}, 
  shape = {kind = 'canoND', 
           object = 
                 {origin = 
                        {x_min+dx/2.0,y_min+height*0.5,z_min+height*0.5}},
                       },
  time_control = {min= {iter= 1}, max = {iter= tmax}, interval = {iter= 10}},
  output = {format = 'ascii'}
},
{
  label = 'probe_middle', 
  folder = 'tracking/',
  variable = {'pressure_phy','velocity_phy'}, 
  shape = {kind = 'canoND', 
           object = 
                 {origin = 
                        {x_min+length*0.5,y_min+height*0.5,z_min+height*0.5}},
                       },
  time_control = {min= {iter= 1}, max = {iter= tmax}, interval = {iter= 10}},
  output = {format = 'ascii'}
},
{
  label = 'probe_outlet', 
  folder = 'tracking/',
  variable = {'pressure_phy','velocity_phy'}, 
  shape = {kind = 'canoND', 
           object = 
                 {origin = 
                        {x_max-dx/2.0,y_min+height*0.5,z_min+height*0.5}},
                      },
  time_control = {min= {iter= 1}, max = {iter= tmax}, interval = {iter= 10}},
  output = {format = 'ascii'}
},
{
  label = 'press_line', 
  folder = 'tracking/',
  variable = {'pressure_phy','press_an', 'diff_pressure'}, 
  shape = {kind = 'canoND', object = {
  origin = {x_min-dx/2.,y_min+height*0.5-offset,z_min+height*0.5-offset},
  vec = {length+2.*dx,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_pressure','press_an',},
--  shape = {kind = 'canoND', object = {
--    origin = {x_min-dx/2.,y_min+height*0.5-offset,z_min+height*0.5-offset},
--    vec = {length+2.*dx,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 = 'wss_spatial_yz',
--  variable = { 'wss_phy', 
--               'wss_an', 
--              'diff_wss'
--             },
--  shape = {kind = 'canoND', object = {
--    origin = {x_min+length*0.5,y_min-offset-dx/2.0,z_min-offset-dx/2.0},
--    vec = {{0.0,0.0,height+2.*dx}, {0.0,height+2.*dx,0.0}},
--    segments = {nHeight+2,nHeight+2} }},
--  time_control = { min = {iter=tmax}, max = {iter=tmax}, interval = {iter=tmax+1} },
--  folder = './tracking/',
--  output = {format = 'asciiSpatial'}
--},
--
--{ label = 'wss_l2norm_yz',
--  variable = {
--    'diff_wss',
--    'wss_an'
--              },
--  shape = {kind = 'canoND', object = {origin ={x_min+length*0.5,y_min-offset-dx/2.0,z_min-offset-dx/2.0},
--                                     vec = {{0.0,0.0,height+2.*dx},
--                                            {0.0,height+2.*dx,0.0}},
--                                     segments = {nHeight+2,nHeight+2} }},
--  time_control = { min = {iter=tmax}, max = {iter=tmax}, interval = {iter=tmax+1} },
--  reduction = {'l2norm', 'l2norm'},
--  folder = './tracking/',
--  output = {format = 'asciiSpatial'}
--},
--{ label = 'vel_spatial_yz',
--  variable = {'vel_mag_phy','velocity_phy'},
--  shape = {kind = 'canoND', object = {origin ={x_min+length*0.5,y_min-dx/2.0,z_min-dx/2.0},
--                                     vec = {{0.0,height+2.*dx,0.0},
--                                            {0.0,0.0,height+2.*dx}},
--                                     segments = {nHeight+2,nHeight+2} }},
--  time_control = { min = {iter=tmax}, max = {iter=tmax}, interval = {iter=tmax+1} },
--  folder = './tracking/',
--  output = {format = 'asciiSpatial'}
--},
--{ label = 'vel_l2norm_yz'}
--  variable = {
--    'diff_vel',
--    'vel_an'
--              },
--  shape = {kind = 'canoND', object = {origin ={x_min+length*0.5,y_min-dx/2.0,z_min-dx/2.0},
--                                     vec = {{0.0,height+2.*dx,0.0},
--                                            {0.0,0.0,height+2.*dx}},
--                                     segments = {nHeight+2,nHeight+2} }},
--  time_control = { min = {iter=tmax}, max = {iter=tmax}, interval = {iter=tmax+1} },
--  reduction = {'l2norm', 'l2norm'},
--  folder = './tracking/',
--  output = {format = 'asciiSpatial',
--},
{ label = 'vel_centerline_x',
  variable = {'vel_mag_phy','velocity_phy','diff_vel'},
  shape = {kind = 'canoND', object = {origin ={x_min-dx/2.0,y_min+height*0.5,z_min+height*0.5},
                                     vec = {length+2.*dx,0.0,0.0},
                                     segments = nLength+2 }},
  time_control = { min = {iter=tmax}, max = {iter=tmax}, interval = {iter=tmax+1} },
  folder = './tracking/',
  output = {format = 'asciiSpatial'}
},
{ label = 'vel_height',
  variable = {'vel_mag_phy','velocity_phy'},
  shape = {kind = 'canoND', object = {origin ={length*0.5,y_min-dx/2.0,0.},
                                     vec = {0.0,height+dx,0.0},
                                     segments = nHeight+2}},
  time_control = { min = {iter=tmax}, max = {iter=tmax}, interval = {iter=tmax+1} },
  folder = './tracking/',
  output = {format = 'asciiSpatial'}
},

}
--print("press_drop = ",press_drop)
--print('press_inlet = ',press_inlet)
--print("u_max for u_inflow = ",u_analy(0.,0.,0.,0.))
--print("u_max_L for u_inflow = ",u_in_L)
--print("Ma_L = ",u_in_L/cs)
--print("wall shear stress at the outside = ",wssA(0.,math.sqrt(2.)*radius,math.sqrt(2.)*radius,0.))
--print('wssA(0,0,0,0) = ',wssA(0.,0.,0.,0.))
