-- Do not alter this file.
-- It is part of the testsuite for regression check
--
require "seeder"

intp         = 'quadratic'
wall_kind    = 'wall_libb'
wall_kindObs = 'wall_linearInterpolation'
inlet_kind   = 'outlet_pab'
outlet_kind  = 'outlet_pab'

scaling = 'diffusive'
scaling = 'acoustic'
relaxation = 'bgk'
stencil    = 'd2q9'
model      = 'lbm_incomp'
Re         = 10.0
initChannel = true
interpolation_method = { method = intp }

tracking_folder = './tracking/'
mesh            = './mesh/'
printRuntimeInfo = false
timing_file = 'timing.res'
simulation_name = 'channelIncli'
control_routine = 'fast'
remove_solid = true

identify = {
  layout = stencil,
  kind   = model,
  relaxation = relaxation,
}

logging = { level = 10 }

-- Reference physical values
u0Phys=1.
viscPhys = u0Phys*heightPhys/Re
csPhys = 300    -- m/s
rho0Phys = 1.
-- Reference LB values
csLB = 1/math.sqrt(3)
cs2=1./3.
u0LB=0.01
rho0LB=1.
omega0=1.7
-- determine the relationships
if scaling == 'acoustic' then
  uLB=u0LB
  dt = uLB/u0Phys*dx
  viscLB = viscPhys*dt/dx/dx
  omega = 1./(3.*viscLB + 0.5)
else
  omega = omega0
  viscLB  = 1./6.*(2./omega - 1.)
  dt = viscLB/viscPhys*dx*dx
  uLB=u0Phys*dt/dx
end
-- Reynolds number
Re = heightPhys*u0Phys/viscPhys
ReLB = heightLB*uLB/viscLB

tEnd = 10
-- tEnd = dt

amplitude = u_in
p0 = rho0LB*cs2*dx*dx/dt/dt

originX = -length*0.25
originY = height*0.
originZ = 00.0
halfwidth = length*0.015
amplitude = 0.00
background = p0
function ic_2Dgauss_pulse(x, y, z)
  return background+amplitude*math.exp(-0.5/(halfwidth^2)*(( x - originX )^2+( y - originY )^2))
end

-- Reference values for the flow state in the 2d stationary channel at laminar
-- flow state
function velArrayRef(x,y,z,t)
  yRotated = -x*math.sin(angle/180*math.pi) + y*math.cos(angle/180*math.pi) 
  uPhys = u0Phys*(1-(2.*yRotated/heightPhys)^2)
  return { uPhys, 0., 0. }
end

function uX(x,y,z,t)
  if walls then
    yRotated = -x*math.sin(angle/180*math.pi) + y*math.cos(angle/180*math.pi) 
    uPhys = u0Phys*(1-(2.*yRotated/heightPhys)^2)
  else
    uPhys = 0.0
  end
  return uPhys
end

function pressureRef(x,y,z,t)
  if walls then
    dp = u0Phys*8.*viscPhys*rho0Phys/heightPhys^2*length
    xRotated = x*math.cos(angle/180*math.pi) + y*math.sin(angle/180*math.pi) 
    return p0 + dp*0.5 - dp/length*xRotated
  else
    dp = 0
    --return ic_2Dgauss_pulse(x,y,z)
    return 0.
  end
end
-- pressureIn=pressureRef(-length/2,0.,0.)
-- pressureOut=pressureRef(length/2,0.,0.)
function Sxx(x,y,z)
  return 0.
  end
function Syy(x,y,z)
  return 0.
  end
function Sxy(x,y,z)
  if walls then
    yRotated = -x*math.sin(angle/180*math.pi) + y*math.cos(angle/180*math.pi) 
    tauxy= -viscPhys*rho0Phys*8./heightPhys^2*u0Phys*yRotated
    S_xy = tauxy/viscPhys/rho0Phys
    return S_xy
  else 
    return 0.
    end
end
function stressRef(x,y,z,t)
  if walls then
    return Sxy(x,y,z)*rho0Phys*viscPhys
  else
    return 0.
  end
end

-- Consistent initial conditions for the channel
if initChannel then
  tRampMin = 1.
  function ic_uX(x,y,z)
    return uX(x,y,z,0)
  end
  function ic_pressure(x,y,z)
    if( model=='lbm') then
      return pressureRef(x,y,z,0) + rho0Phys*cs2*dx/dt
    else
      return pressureRef(x,y,z,0)
    end
  end
  function ic_Sxx(x,y,z)
    return Sxx(x,y,z)
  end
  function ic_Syy(x,y,z)
    return Syy(x,y,z)
  end
  function ic_Sxy(x,y,z)
    return Sxy(x,y,z)
  end
else
  tRampMin = 0.
  function ic_uX(x,y,z)
    return 0.0
  end
  function ic_pressure(x,y,z)
    return ic_2Dgauss_pulse(x,y,z)
  end
  function ic_Sxx(x,y,z)
    return 0.
  end
  function ic_Syy(x,y,z)
    return 0.
  end
  function ic_Sxy(x,y,z)
    return 0.
  end
end

variable = {
  -- analytical solutions
  { name='pressureRef', ncomponents=1, vartype = 'st_fun', st_fun = pressureRef },
  { name='pressureIn',  ncomponents=1, vartype = 'st_fun', st_fun = pressureRef(-length/2,0.,0.) },
  { name='pressureOut', ncomponents=1, vartype = 'st_fun', st_fun = pressureRef( length/2,0.,0.) },
  -- { name='pressureIn',  ncomponents=1, vartype = 'st_fun', st_fun = 0 },
  -- { name='pressureOut', ncomponents=1, vartype = 'st_fun', st_fun = 0 },
  { name='velMagRef',   ncomponents=1, vartype = 'st_fun', st_fun = uX },
  { name='velArrayRef', ncomponents=3, vartype = 'st_fun', st_fun = velArrayRef },
  { name='stressRef',   ncomponents=1, vartype = 'st_fun', st_fun = stressRef },
  -- error
  { name='error_v', ncomponents=1, vartype = 'operation',
    operation = {
      kind = 'difference',
      input_varname = {'vel_mag_phy', 'velMagRef'}
    },
  },
}
init_allElems = false
-- Time step settigs
interval = tEnd/10.
tRamping = tEnd/10.
sim_control = {
  time_control = { max = tEnd, interval = interval }
}
physics = { dt = dt, rho0 = rho0Phys }

fluid = {
  omega = omega
}

-- Initial condition in physical units
initial_condition = {
  pressure  = ic_pressure,
  velocityX = ic_uX,
  velocityY = 0.0,
  velocityZ = 0.0,
  Sxx = ic_Sxx,
  Syy = ic_Syy,
  Sxy = ic_Sxy,
}
-- Boundary conditions
boundary_condition = {
  { label = 'west', kind = inlet_kind,  pressure = 'pressureIn',  },
  { label = 'east', kind = outlet_kind, pressure = 'pressureOut', },
  { label = 'north', kind = wall_kind },
  { label = 'south', kind = wall_kind },
}


if usePeriodic == false then
  table.insert( boundary_condition,
{ label = 'frontback',
   kind = 'wall' } )
end

if useObstacle == true then
  table.insert( boundary_condition,
{ label = 'sphere',
   kind = wall_kindObs } )
end


-- Tracking
tracking = {
  {
    label = 'probePressure_l'..level,
    variable = {
     'pressure_phy',
     'velocity_phy',
     'shear_stress_phy',
    },
    shape = {
      kind = 'canoND',
      object = {
         origin = {  -math.sin(angle/180*math.pi)*length*0.,
                     -math.cos(angle/180*math.pi)*height*0.,
                    0.5*dx},
       },
    },
    time_control = {
      useIterations = true, min = 1, max = math.floor(tEnd/dt), interval = 1
    },
    output={format = 'ascii'}, folder = tracking_folder,
  },

 {
    -- tracking object for getting the error in terms of the l2norm
    label = 'global_l'..level,
    variable = {
       'pressure_phy',
       'velocity_phy',
       'shear_stress_phy',
       -- 'pdf',
       'treeid',
     },
    shape = { kind = 'all', },
    time_control = {min = tEnd, max = tEnd, interval = tEnd},
    output={format = 'vtk'},
    folder = tracking_folder,
 },
 {
  -- tracking object for getting the error in terms of the l2norm
  label = 'dpdx'..level,
  variable = {
    'pressure_phy', 'pressureRef'
  },
  shape = {
    kind = 'canoND',
    object = {
      origin ={ -length*0.25,-height*0.,0.5*dx},
      vec={
        {length*0.6, 0., 0.},
        {0.,0.,0.}
      },
      segments = {2*nElemsMax, nElemsMax/2},
    },
  },
  time_control = {min = tEnd, max = tEnd, interval = interval},
  output={format = 'asciiSpatial'},
  folder = tracking_folder
 },

 {
  -- tracking shear stress
  label = 'shear_stress_'..level,
  variable = { 
    'shear_stress_phy', 'stressRef',
  },
  shape = {
    kind = 'canoND',
    object = {
      origin ={ -length*0.,-height*0.5,0.5*dx}, 
      vec={{length*0., height, 0.}, {0.,0.,0.}}, 
      segments = {2*nElemsMax, 1}
    },
  },
  time_control = {min = tEnd, max = tEnd, interval = interval},
  output={format = 'asciiSpatial'},
  folder = tracking_folder
 },

 {
  -- tracking object for getting the error in terms of the l2norm
  label = 'velErr'..level,
  variable = {
    'error_v',
  },
  reduction = {'l2norm'},
  shape = {
    kind = 'canoND',
    object = {
      origin = {  -math.sin(angle/180*math.pi)*length*0.,
                  -math.cos(angle/180*math.pi)*height*0.5,
                 0.5*dx},
      vec = {
        {-height*math.sin(angle/180*math.pi), height*math.cos(angle/180*math.pi), 0.},
        {-height*math.sin(angle/180*math.pi), height*math.cos(angle/180*math.pi), 0.}
      },
      segments = {2*nElemsMax, nElemsMax},
    },
  },
  time_control = {min = tEnd, max = tEnd, interval = dt},
  output={format = 'asciiSpatial'},
  folder = tracking_folder
 },

 {
  -- track velocity across channel
  label = 'vel'..level,
  variable = {
    'error_v',
    'vel_mag_phy',
    'velMagRef',
  },
  shape = {
    kind = 'canoND', 
    object = {
      origin = {  -math.sin(angle/180*math.pi)*length*0.,
                  -math.cos(angle/180*math.pi)*height*0.5,
                 0.5*dx}, 
      vec = {
        {-height*math.sin(angle/180*math.pi), height*math.cos(angle/180*math.pi), 0.},
        {-height*math.sin(angle/180*math.pi), height*math.cos(angle/180*math.pi), 0.}
      }, 
      segments = {2*nElemsMax, nElemsMax} 
    } 
  },
  time_control = {min = tEnd, max = tEnd, interval = interval},
  output={format = 'asciiSpatial'},
  folder = tracking_folder
 },

} -- tracking

estart = {
  write = './restart/',
  time_control = {min = tEnd, max = tEnd, interval = tEnd }
}

if verbose then
  print('           omega   = '..omega)
  print('          viscLB   = '..viscLB)
  print('        heightLB   = '..heightLB)
  print('             uLB   = '..uLB)
  print('Reynolds number Re = '..ReLB)
  print('          viscPh   = '..viscPhys)
  print('        heightPh   = '..heightPhys)
  print('        lengthPh   = '..length)
  print('             uPh   = '..u0Phys)
  print('Reynolds number Re = '..Re)
  -- print( '  pIn '.. pressureIn )
  -- print( '  pOut'.. pressureOut)
end

ebug = {
  verbosity = 100,
  logging = {
    level = 10,
    filename = 'dbg',
    root_only = true,
    real_form = 'EN12.3',
    int_form  = 'I0',
  },
  debugMode        = true,
  debugFiles       = true,
  traceMemory      = false,
  dumpDependencies = false,
  dumpTreeIDs      = false,
  dumpAuxLists     = false,
  checkEachStep    = false,
  dumpBoundaries   = false,
}

