-- Musubi configuration file. 
require 'common'

--flow parameters
rho0_L = 1.0
cs = math.sqrt(1.0/3.0)

MachNr = 0.01
u_fric_L = MachNr*cs
dt = u_fric_L*dx/u_fric_p
nu_L = nu_p*dt/dx^2
omega = 1.0/(3.0*nu_L+0.5)

-- reference pressure
ref_press = rho_p*dx^2/dt^2

bc_north = 'symmetry'
bc_south = 'wall'

function u_testgrad(x,y,z) 
  return (0.5/(rho_p*nu_p))*dpdx*y*(height-y)
end



-- Paper: Simons, T., & Pletcher, R. (2013). Large eddy simulation of turbulent 
-- flows using unstructured grids, (October). https://doi.org/10.2514/6.1998-3314
-- magnitude of perturbation
mag_pert = 3.65e-5
PI = math.pi
-- perturbation in x-velocity
function u_x_perturbation(x,y,z)
  -- normalized distance
  x_star = x/length
  y_star = y/height
  z_star = z/width 
  term1 = math.cos(2*PI*x_star)*math.sin(2*PI*y_star)
  term2 = 0.5*math.cos(4*PI*x_star)*math.sin(2*PI*y_star)
  term3 = math.cos(2*PI*x_star)*math.sin(4*PI*y_star)
  return mag_pert*length*math.sin(PI*z_star)*(term1+term2+term3)
end 

-- perturbation in y-velocity
function u_y_perturbation(x,y,z)
  -- normalized distance
  x_star = x/length
  y_star = y/height
  z_star = z/width 
  term1 = 0.5*math.sin(2*PI*x_star)*math.cos(2*PI*y_star)
  term2 = 0.5*math.sin(4*PI*x_star)*math.cos(2*PI*y_star)
  term3 = 0.25*math.sin(2*PI*x_star)*math.cos(4*PI*y_star)
  return -mag_pert*height*math.sin(PI*z_star)*(term1+term2+term3)
end 

-- perturbation in z-velocity
function u_z_perturbation(x,y,z)
  -- normalized distance
  x_star = x/length
  y_star = y/height
  z_star = z/width 
  term1 = math.sin(2*PI*x_star)*math.sin(2*PI*y_star)
  term2 = math.sin(4*PI*x_star)*math.sin(2*PI*y_star)
  term3 = math.sin(2*PI*x_star)*math.sin(4*PI*y_star)
  return -mag_pert*(1.0+math.cos(PI*z_star))*(term1+term2+term3)
end 

--karman constant
k = 0.41
function u_x_initial(x,y,z)
  if (y>halfHeight) then
    y_plus = u_fric_p*(height-y)/nu_p
  else
    y_plus = u_fric_p*y/nu_p
  end
  y_plus_div11 = y_plus/11.0
  logterm = math.log(1.0+0.4*y_plus)/k
  expterm = math.exp(-y_plus_div11)
  U_R = logterm + 7.8*(1.0 - expterm -y_plus_div11 * expterm)
  u_x = u_fric_p * U_R
  --return u_x + (math.random()*2.0-1.0)*u_fric_p*0.05
  return u_x + u_x_perturbation(x,y,z) 
end

function u_initial_rand(x,y,z)
  return (math.random()*2-1)*u_fric_p*0.05
end

-- Simulation name
simulation_name = 'Channel'
mesh = 'mesh/' -- Mesh information
printRuntimeInfo = false
scaling = 'acoustic'
interpolation_method = 'quadratic'
NOdebug = {logging = {level=5, filename='dbg', root_only=false}}
logging = {level=5}
-- Time step settigs
tmax_p = 1*T_c --s
tmax =  tmax_p/dt    -- total iteration number
interval_p = 1 -- s
interval = interval_p/dt
sim_control = {
  time_control = { 
    max = {sim=tmax_p,iter=1}, 
    interval = {sim=tmax_p/100}
   },
  abort_criteria = {
    stop_file = 'stop',
  }
}

restart = {
  NOread = 'restart/Channel_lastHeader.lua',
  write = 'restart/',
}

fluid = { 
  kine_shear_viscosity = nu_p, -- old implementation
  kinematic_viscosity = nu_p,
  bulk_viscosity = nu_p,
  rho0 = rho_p 
}


-- Initial condition 
initial_condition = { 
  pressure = ref_press, 
  velocityX = u_x_initial,
  velocityY = u_y_perturbation,
  velocityZ = u_z_perturbation 
}

identify = {
  label='turb',
  layout='d3q19', 
  relaxation='bgk_les_new',
  --relaxation = 'bgk_les',
  kind = 'lbm'
}

c_s = 0.17
turbulence = {
  model = 'vreman',
  c_s = c_s,
  c_v = 0.5,--math.sqrt(2.5)*c_s,
  c_w = 0.5
}

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

-- Boundary conditions
boundary_condition = {  
  {   
    label = 'north', 
    kind = bc_north,
  },
  {
    label = 'south', 
    kind = bc_south,
  },
}

source = {
  force = {dpdx,0.0,0.0}
}

NOvariable = {
  { name='press_avg', 
    ncomponents=1, 
    vartype = 'operation', 
    operation = {
      kind = 'reduction_transient',
      input_varname = {'pressure_phy'},
      reduction_transient = {
        kind = 'average',
        nrecord = 14*T_c
      }
    }
  },
  { name='velocity_avg', 
    ncomponents=3, 
    vartype = 'operation', 
    operation = {
      kind = 'reduction_transient',
      input_varname = {'velocity_phy'},
      reduction_transient = {
        kind = 'average',
        nrecord = 14*T_c
      }
    }
  },
}

-- Tracking              
tracking = {
  {
    label = 'vtk', 
    folder = 'tracking/',
    NOvariable = {'pressure_phy','velocity_phy', 'press_avg','velocity_avg'}, 
    variable = {'pressure_phy','velocity_phy','turb_viscosity_phy','kine_viscosity_phy'}, 
    shape = {kind = 'all'},
    time_control = {min= 0, max = tmax_p, interval = tmax_p/50},
    output = {format = 'vtk'}
  },
  {
    label = 'probe', 
    folder = 'tracking/',
    NOvariable = {'pressure_phy','velocity_phy', 'press_avg','velocity_avg'}, 
    variable = {'pressure_phy','velocity_phy'}, 
    shape = {
      kind = 'canoND',
      object = {origin = {0.0,halfHeight-dx/2.0,0.0},
                vec = { {length,0.0,0.0},
                        {0.0,0.0,width} }
      }
    },
    reduction = {'average','average'},
    time_control = {min= 0, max = tmax_p, interval = 100*dt},
    output = {format = 'ascii'}
  },
}
