require "seeder"

intp = 'linear'
intp = 'compact'

wall_kind = 'wall_linearInterpolation'
wall_kind = 'wall'

inlet_kind  = 'outlet_expol'
outlet_kind = 'outlet_expol'

scaling = 'diffusive'
scaling = 'acoustic'

relaxation = 'mrt'
stencil    = 'd2q9'

model = 'lbm_incomp'
model = 'lbm'

tEnd = 10
Re   = 60

initChannel = false
testIntp = true
control_routine = true
track = false
restart = false
recheck = true
global  = true

interpolation_method = {
  method = intp,
}

--require "PATH"
tracking_folder = './tracking/'--"TRACKING_FOLDER"
output_folder = './output/'--"OUTPUT_FOLDER"

-- Simulation name
mesh = './mesh/' -- Mesh information
simulation_name = simName

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

-- Reference physical values
u0Phys=1.
viscPhys = u0Phys*heightPhys/Re
csPhys   = 300    -- m/s
rho0Phys = 1.
-- Reference LB values
csLB    = 1.0 / math.sqrt(3.0)
cs2     = 1.0 / 3.0
u0LB    = 0.05
rho0LB  = 1.0
rho0LB0 = 0.0
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 = (2/omega - 1.) / 6
  dt     = viscLB/viscPhys*dx*dx
  uLB    = u0Phys*dt/dx
end

-- Reynolds number
Re = heightPhys*u0Phys/viscPhys
ReLB = heightLB*uLB/viscLB
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)
end


amplitude = u_in
if model == 'lbm_incomp' then
  p0 = 0.0
else
  p0 = rho0LB*cs2*dx*dx/dt/dt
  p0 = 0
end

originX = -length*0.25
originY = height*0.
originZ = 0.0
halfwidth = length*0.015
amplitude = 0.0
background = p0

function ic_2Dgauss_pulse(x, y, z, t)
  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 uX(x,y,z, t)
    if walls then
      uPhys = u0Phys*(1-(2.*y/heightPhys)^2)
    else
      uPhys = 0. --u0Phys
    end
    return uPhys
end

function pressureRef(x,y,z, t)
  if walls and periodicX == false then
    dp = u0Phys*8.*viscPhys*rho0Phys/heightPhys^2*length
    if useObstacle or ( not usePeriodic or not case2d ) then
      dp = 10.*dp
    end
    return p0 + dp*0.5 - dp/length*x
  else
    dp = 0
    return 0.
  end
end
pressureIn  = pressureRef(-length/2,0.,0.)
pressureOut = pressureRef(length/2,0.,0.)
function Sxx(x,y,z, t)
  return 0.
  end
function Syy(x,y,z,t)
  return 0.
  end
function Sxy(x,y,z,t)
  if walls then
    tauxy= -viscPhys*rho0Phys*8./heightPhys^2*u0Phys*y
    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,t)*rho0Phys*viscPhys
  else
    return 0.
  end
  end

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

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_ramping = { 
  --  predefined='smooth', 
  -- min_factor = 0.7, 
  -- max_factor=1, 
  -- from_time=0, 
  -- to_time=tRamping}, 
       omega = omega, 
       rho0 = rho0Phys,
  }

-- 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 = {  
 }
if walls then
  table.insert( boundary_condition, 
{ label = 'north', 
   kind = 'wall' })
  table.insert( boundary_condition, 
{ label = 'south', 
   kind = 'wall' })
end 


if usePeriodic ==false then
  table.insert( boundary_condition, 
{ label = 'frontback', 
   kind = 'wall' } )
    if periodicX == false then
      table.insert( boundary_condition,
    { label = 'west', 
      kind = inlet_kind,
      pressure = pressureIn } )
    end
else
  if periodicX == false then
    table.insert( boundary_condition,
  { label = 'west', 
    kind = inlet_kind,
    pressure = pressureIn, 
  } )
  end
end

if periodicX == false then
  table.insert( boundary_condition, 
{ label = 'east', 
   kind = outlet_kind, pressure = pressureOut })
end

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

variable = {
  -- reference results
  { name='pressureRef', ncomponents=1, vartype = 'st_fun', st_fun=pressureRef },
  { name='stressRef',   ncomponents=6, vartype = 'st_fun', st_fun=stressRef },

  -- error
  { name='error_p',  ncomponents=1, vartype = 'operation',
    operation = {
      kind = 'difference',
      input_varname = {'pressure_phy', 'pressureRef'}
    },
  },
}

if track then
-- Tracking
tracking = {

  -- flow rate measurement at the inlet
  { label = 'flowR', 
    variable = {'velocity_phy'}, 
    reduction = {'average'},
    shape = {
      kind = 'canoND', 
      object = {
        origin ={
          -length*0.5+3*dx,
          -height/2, -height/2
        },
        vec = {
          {0., height, 0.},
          {0., 0., height}
        },
        segments = {nElemsMax/2, nElemsMax/2}
      }
    },
    output={format='ascii'}, folder=tracking_folder, 
    time_control = { interval = interval, min = 0, max = tEnd }
  },

  -- probe at point: 0 0 0
  {
    label = 'probePressure_l'..level, 
    variable = {
      'pressure_phy',
      'velocity_phy',
      'shear_stress_phy',
    },
    shape = {
      kind = 'canoND',
      object = {
        origin ={0.0,0.,0.}
      }
    },
    time_control = { min = 1, max = tEnd, interval = 1},
    output={format = 'ascii'}, folder = tracking_folder,
  },
-- {
--  -- tracking object for getting the error in terms of the l2norm
--  label = 'hvsXZ_l'..level, 
--  variable = { {'pressure_phy'},
--     {'velocity_phy'}, {'shearstress_phy'}},
--  shape = {kind = 'canoND', object = {origin ={ -length*0.5,0.,-height*0.5}, 
--            vec={{length, 0., 0.}, {0.,0.,height}}, 
--            segments = {2*nElemsMax, nElemsMax/2} } },
--  time_control = {min = 0., max = tEnd, interval = interval},
--  output={format = 'harvester'},
--  folder = tracking_folder
-- }, 
  {
  -- tracking object for getting the error in terms of the l2norm
  label = 'hvsXY_l'..level, 
  variable = { 
     'pressure_phy', 
     'velocity_phy',
     'shear_stress_phy',
   },
  shape = {
    kind = 'canoND',
    object = {
      origin ={ -length*0.5,-height*0.5,0.5*dx}, 
      vec={{length, 0., 0.}, {0.,height,0.}}, 
      segments = {2*nElemsMax, nElemsMax/2} } },
  time_control = {min = 0., max = tEnd, interval = interval},
  output={format = 'harvester'},
  folder = tracking_folder
 },

 -- tracking pressure along channel
 {
  label = 'dpdx'..level, 
  variable = { 
     'pressure_phy',
     'pressureRef',
  },
  shape = {
    kind = 'canoND',
    object = {
      origin ={ -length*0.25, 0, dx/2}, 
      vec = {
        {length*0.6, 0., 0.},
      },
      segments = {2*nElemsMax }
    }
  },
  time_control = {min = tEnd, max = tEnd, interval = interval},
  output={format = 'asciiSpatial'},
  folder = tracking_folder,
 },

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

 -- tracking pressure error along channel
 {
   label = 'Errdpdx'..level, 
   variable = { 
     'error_p',
   },
   reduction = {'l2norm'},
   shape = {
     kind = 'canoND',
     object = {
       origin ={ -length*0.25,0,dx/2}, 
       vec={{length*0.6, 0., 0.} }, 
       segments = {2*nElemsMax }
     }
   },
   time_control = {min = tEnd, max = tEnd, interval = interval},
   output={format = 'ascii'},
   folder = tracking_folder,
 },
} -- tracking table
else
  tracking = {}
end

if recheck then
  table.insert( tracking,
  {
    label = 'reference',
    variable = {'pressure_phy'},
    shape = {
      kind = 'canoND',
      object = {
        origin = { 0.0 , 0.0, 0.0 },
      },
    },
    time_control = {min = tEnd, max = tEnd, interval=tEnd},
    output={format = 'ascii'}, folder = tracking_folder,
   }
   )
end

if global then
  table.insert( tracking,
  {
    label = 'global',
    variable = {
      'pressure_phy',
      'velocity_phy',
      -- 'strain_rate_phy',
    },
    shape = {
      kind = 'global',
    },
    time_control = {min = tEnd, max = tEnd, interval=tEnd/10},
    output={format = 'vtk'}, folder = output_folder,
   }
   )
end

if restart then
restart = {
  write = './restart/',
  time_control = {min = tEnd, max = tEnd, interval=tEnd},
}
end
